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A.  SCIENTIFIC  OBJECTIVES:  The  objective  of  this  research  program  is  to  devise 
innovative  joint  time-frequency  (JTF)  processing  concepts  for  radar  image  enhancement 
and  physics-based  feature  extraction.  In  particular,  we  investigate  how  JTF  techniques 
can  be  utilized  to  enhance  synthetic  aperture  radar  (SAR)  and  inverse  synthetic  aperture 
radar  (ISAR)  imageries  by  removing  artifacts  due  to  uncompensated  target  motion, 
conplex  target  scattering  physics,  articulating  target  components  and  clutter  and 
propagation  effects.  Furthermore,  we  set  out  to  re-interpret  the  extracted  artifacts  in  a 
more  meaningful  feature  space  so  that  they  can  be  utilized  to  enhance  the  performance  of 
target  identification  algorithms.  This  research  is  leveraged  against  our  previous  JTF  work 
under  the  Joint  Services  Electronics  Program,  as  well  as  our  state-of-the-art  radar 
signature  simulation  capabilities.  The  JTF-based  radar  image  processing  took  developed 
under  this  program  will  be  disseminated  to  radar  researchers  in  the  US  Navy. 

B.  SUMMARY  OF  RESULTS  AND  SIGNEPICANT  ACCOMPLISHMENTS: 

During  this  prograno,  we  have  followed  two  lines  of  research:  (i)  processing  real  radar 
data  to  identify  needed  areas  of  research  and  to  test  our  JTF  algorithms,  (ii)  developing 
new  JTF  algorithms  to  address  the  problem  areas  identified  in  (i). 

Our  accomplishments  along  the  first  hne  of  research  include: 


1 . 1  Application  of  JTF  motion  compensation  algorithm  to  NATO  TIRA  radar  data. 

1 .2  Processing  of  NATO  MERIC  radar  data. 

1 .3  Processing  of  Navy  SCATR  data. 

Along  the  second  line  of  research,  we  have  focused  on  two  key  research  areas 
identified  from  processing  of  measurement  data.  The  first  area  focuses  on  the  imaging  of 
targets  with  movable  components.  We  have  developed  a  number  of  JTF  algorithms  to 
extract  features  from  targets  having  moving  components: 

2. 1  Removal  of  image  artifacts  due  to  helicopter  rotor  blades. 

2.2  Removal  of  image  artifacts  due  to  jet  engine  modulation. 

2.3  Extraction  of  microDoppler  features  from  ISAR  data  using  adaptive  chirplet 
representation. 

The  second  area  focuses  on  targets  having  highly  chaotic  motions.  We  have 
addressed  a  number  of  issues  in  image  formation  in  the  presence  of  chaotic  target 
motions: 

3. 1 .  Three-dimensional  motion  detection  using  JTF  algorithm. 

3.2.  3D  ISAR  image  reconstruction  of  a  target  with  motion  data. 

3.3.  Application  of  genetic  algorithms  to  ISAR  imaging. 

Finally,  we  have  carried  out  some  exploratory  research  on  two  topics  related  to  radar 
imaging  and  JTF  analysis: 

4. 1 .  Clutter  reduction  for  SAR  images  using  adaptive  wavelet  packet  transform. 

4.2.  Inverse  scattering  using  genetic  algorithms. 

The  detailed  descriptions  of  our  accomplishments  are  discussed  below. 

1.1.  Application  of  JTF  Motion  Compensation  Algorithm  to  NATO  TIRA  Data. 

We  have  devoted  considerable  efforts  to  process  the  NATO  TIRA  radar  data  taken  in 
Germany  in  November  1997.  The  objectives  of  this  effort  are  to  test  our  previously 
developed  adaptive  joint  time-frequency  (AJTF)  algorithm  [1]  for  ISAR  motion 
compensation  and  to  identify  needed  areas  of  research  in  ISAR-based  target  recognition. 
Our  algorithm  uses  a  search  and  projection  technique  in  the  joint  (dwell  time)-(Doppler 
frequency)  plane  to  select  and  track  the  prominent  point  scatterers.  The  higher-order 
translation  and  rotation  motions  are  then  extracted  and  compensated  for  in  the  data  to 


form  a  focused  image  of  the  target.  The  motion  compensated  images  have  been 
compared  against  the  reference  images  generated  by  using  the  motion  information 
available  in  the  instrumented  ARDS  data.  Furthermore,  comparison  has  also  been  made 
with  the  simulated  ISAR  images  of  the  air  target  from  the  radar  signature  prediction  code 
Xpatch.  A  set  of  representative  images  is  shown  below.  Figs.  1.1.1(a),  (b)  and  (c)  are  the 
images  for  a  target  at  azimuth=56°  (0°  being  nose-on)  generated  using,  respectively, 
AJTF  motion  compensation,  ARDS  sensor  information  and  Xpatch  simulation.  The 
dynamic  range  of  the  displayed  images  is  55  dB.  The  look  angle  information  was 
deciphered  from  the  ARDS  data.  By  comparing  Figs.  1.1.1(a)  and  1.1.1(b),  we  observe 
that  the  motion  compensated  image  and  the  ARDS-derived  reference  image  appear  to  be 
in  good  agreement.  The  Xpatch  image  is  more  focused  and  does  not  exhibit  the  diffused 
characteristics  of  the  measurement  data.  A  more  detailed  CAD  model  should  improve 
the  quality  of  the  predicted  image. 

While  the  comparison  among  the  images  is  quite  good,  we  have  identified  several 
needed  areas  of  research  in  ISAR-based  target  recognition.  Fig.  1.1.2  shows  both  the 
correlation  coefficient  between  the  JTF  and  tmth  images  and  that  between  the  JTF  and 
Xpatch  images  versus  azimuth  look  angle.  From  the  two  curves,  we  see  that  the  JTF 
images  agree  very  well  with  the  truth  images,  indicating  that  blind  motion  compensation 
is  a  very  feasible  method  for  processing  real-world  radar  data.  The  correlation  coefficient 
between  the  JTF  and  Xpatch  images  is  slightly  lower  than  that  between  the  measured 
images.  In  particular,  two  problem  regions  can  be  seen  from  this  plot.  First,  in  the  region 
near  nose-on  (180  degrees  in  azimuth),  the  correlation  coefficient  is  significantly  lower. 
The  reason  is  due  to  the  strong  jet  engine  modulation  (JEM)  lines  in  the  measured  data. 
This  problem  is  further  discussed  in  2.2  and  an  algorithm  to  remove  JEM  lines  is 
proposed.  Second,  at  some  angles  around  the  broadside  region  (90  degrees  in  azimuth), 
the  correlation  coefficient  is  also  low.  In  this  case,  both  the  associated  JTF  image  and  the 
truth  image  are  both  found  to  be  of  low  quality.  After  further  investigation,  it  is  found 
that  the  image  blurring  is  due  to  the  non-steadiness  of  the  imaging  plane.  This  problem  is 
further  addressed  in  3.1. 


1.2.  Processing  of  NATO  MERIC  Radar  Data.  We  have  also  processed  the  MERIC 
radar  data  made  available  to  us  through  the  Naval  Research  Laboratory.  Blind  motion 
compensation  is  more  challenging  in  this  case,  as  the  flight  path  of  the  aircraft  was  very 
close  to  the  radar  and  the  instrument  data  showed  more  variations  in  azimuth  and 
elevation  angles  during  the  flight  of  the  aircraft.  The  task  is  further  complicated  by  the 
lack  of  sufficient  instrument  data  (some  of  the  ARDS  pod  data  had  been  heavily  down- 
sampled  and  written  into  Microsoft  Excel  form  in  the  MERIC  data  set).  In  addition, 
some  of  the  radar  data  contain  very  strong  colored  noise  in  the  background  and 
considerable  effort  was  spent  on  denoising  the  data.  Despite  these  difficulties,  the  final 
JTF  image  quality  we  have  generated  is  fairly  good.  Fig.  1.2.1  shows  the  angular  motion 
parameters  of  an  aircraft  for  which  the  radar  data  are  available.  The  aircraft  undergoes  an 
azimuth  rotation  of  about  25  degrees  during  the  17.5-sec  collection  period.  While  there 
are  35,000  pulses  of  radar  range  profiles  during  this  duration,  only  18  samples  of 
instrument  data  are  available  in  the  same  duration.  Consequently,  no  ground  truthing  was 
possible  due  to  the  lack  of  sufficient  instrument  data.  To  validate  our  results,  we  have 
generated  the  Xpatch  prediction  using  a  CAD  model  purchased  from  viewpoint.com.  It 
is  however  not  the  exact  model  for  the  aircraft  used  in  the  measurement.  The  model  was 
first  converted  into  an  Xpatch-compatible  format.  It  was  then  edited  to  remove  the 
landing  gear  and  bomb  loads  under  the  wings,  as  well  as  to  seal  up  the  open  seams  in  the 
model.  Some  comparison  results  are  presented  in  Fig.  1.2.2.  Figs.  1.2.2(a)  and  1.2.2(b) 
show  respectively  the  JTF  motion  compensated  image  and  the  Xpatch  predicted  image 
near  point  P  of  Fig.  1.2.1  in  the  flight  path.  The  agreement  between  the  MERIC  image 
and  the  Xpatch  simulated  image  is  fairly  good,  considering  the  uncertainties  in  the  CAD 
model.  Figs.  1.2.2(c)  and  1.2.2(d)  show  respectively  the  JTF  motion  compensated  image 
and  the  Xpatch  predicted  image  near  point  Q  of  Fig.  1.2.1  in  the  flight  path.  It  is  clear 
that  the  agreement  is  quite  poor  in  this  case.  The  blind  image  formation  is  quite 
challenging  during  this  particular  portion  of  the  flight.  Furthermore,  we  are  not  certain  of 
the  exact  imaging  plane  for  the  simulation  due  to  the  sparseness  of  instrument  data.  An 
interesting  phenomenon  was  observed  from  the  rnotion  compensated  image  of  a  second 
aircraft.  The  formed  images  show  Doppler  smearing  near  the  nose  region.  After  some 
research  on  this  aircraft,  we  have  determined  that  it  is  very  likely  due  to  movement  from 
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a  mechanically  scanned  antenna  in  the  nose  cone  of  the  aircraft.  This  topic  is 
investigated  in  more  detail  in  2.3. 
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Fig.  1.2.1.  Azimuth  versus  elevation  angles  of  the 

aircraft  with  respect  to  the  radar  during  the 
flight  path. 
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Fig.  1.2.2.  Comparison  of  the  MERIC  image  to  Xpatch  simulation,  (a)  JTF 
motion  compensated  image  at  point  P.  (b)  Xpatch  simulated 
image  at  point  P.  (c)  JTF  motion  compensated  image  at  point  Q. 
(d)  Xpatch  simulated  image  at  point  Q. 


1.3.  Processing  of  Navy  SCATR  Data.  We  have  also  processed  the  small  craft 
(SCATR)  radar  data  made  available  to  us  through  the  Naval  Research  Laboratory.  In 
contrast  to  the  air  data,  the  motion  of  the  small  boat  is  strongly  dictated  by  the  motion  of 
the  ocean  surface.  As  a  result,  the  radar  image  shows  rapid  variation  in  the  imaging  plane 
and  the  interpretation  of  the  images  becomes  more  challenging.  Fig.  1.3.1(a)  shows  the 
model  of  the  boat.  Fig.  1.3.1(b)  shows  the  ISAR  image  formed  using  the  JTF  motion 
compensation  algorithm  during  a  particular  dwell  duration.  The  quality  of  the  image  is 
fair.  Fig.  1.3.1(c)  shows  the  pose  of  the  target  with  respect  to  the  radar  generated  using 
the  data  recorded  by  the  GPS  sensors  carried  on  the  boat.  Because  the  GPS  data  was  not 
sufficient  to  fully  determine  the  pose,  a  number  of  assumptions  were  made  to  generate 
the  pose  information.  As  a  result,  we  were  not  fully  confident  that  the  radar  data  and  the 
motion  data  have  good  correspondence.  Nevertheless,  we  proceeded  to  interpret  the 
motion  data.  It  is  clearly  seen  from  Fig.  1.3.1(c)  that  the  rotational  motion  of  the  target 
over  a  3°  by  6°  AZ-EL  angular  window  follows  a  rather  chaotic  path  and  is  not  well 
confined  to  a  2D  plane.  Fig.  1.3.1(d)  shows  the  projection  of  the  target  model  into  the 
corresponding  imaging  plane.  As  we  can  see,  during  this  particular  (labeled  in  red) 
imaging  interval,  the  image  correspond  approximately  to  the  sideview  of  the  boat. 
However,  during  the  subsequent  collection  intervals,  the  image  plane  changes  quite 
rapidly  to  other  views  of  the  target.  This  issue  is  similar  to  the  problem  we  have 
encountered  in  air  targets.  However,  it  occurs  much  more  frequently  and  the  effect  is 
much  more  dramatic  for  small  boats.  Moreover,  it  becomes  increasingly  difficult  to  find 
a  sufficiently  long  imaging  interval  even  to  form  a  decent  image.  The  issue  of  finding 
good  imaging  intervals  will  be  further  discussed  in  3.1. 

2.1.  Removal  of  Image  Artifacts  due  to  Helicopter  Rotor  Blades.  It  is  well  known 
that  when  rotating  components  exist  on  a  target  such  as  gimbaled  antennas  or  propeller 
blades,  image  artifacts  are  introduced  in  the  Doppler  dimension  of  the  ISAR  image  [2,  3]. 
These  smeared  features  oftentimes  overshadow  the  target  geometrical  features  and  hinder 
the  proper  interpretation  of  the  ISAR  image.  We  have  developed  a  technique  to  remove 
such  Doppler  smear  and  produce  a  clear  ISAR  image  of  the  target  based  on  adaptive  joint 
time-frequency  processing  [4,  5].  The  technique  entails  adaptively  searching  for  the 


Fig.  1.3.1.  Image  from  the  SCATR  data,  (a)  Target  CAD  model,  (b)  Image 
formed  using  the  AJTF  algorithm,  (c)  Target  pose  data  determined 
from  GPS  data,  (d)  Target  model  projected  into  the  corresponding 
imaging  plane. 


linear  chirp  bases  that  best  represent  the  time-frequency  behavior  of  the  signal.  This  is 
accomplished  by  projecting  the  signal  onto  all  possible  chirp  bases  and  finding  the  one 
with  the  maximum  projection  value.  After  the  optimal  basis  is  found,  the  signal 
component  associated  with  this  basis  is  subtracted  from  the  original  signal.  By  iterating 
this  search  procedure,  the  signal  can  be  fully  parameterized  with  a  set  of  chirp  basis 
functions.  Since  the  Doppler  frequency  due  to  the  rotating  component  is  both  larger  and 
more  rapidly  varying  (in  dwell  time)  than  that  from  the  target  body,  the  signal 
components  due  to  the  fast  rotating  part  are  associated  with  those  chirp  bases  having 
large  displacement  and  slope  parameters.  On  the  other  hand,  the  signal  components 
due  to  the  target  body  motion  are  represented  by  those  chirp  bases  with  relatively  small 


(a)  ISAR  image  of  helicopter 


(b)  ISAR  image  after  Doppler  smearing 
due  to  rotating  blades  is  removed 


Fig.  2.1.1.  Adaptive  joint  time- 
frequency  processing  of  ISAR 
imagery  from  a  helicopter. 
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displacement  and  slope  parameters.  By  sorting  these  chirp  bases  according  to  their  slopes 
and  displacements,  the  scattering  due  to  the  fast  rotating  part  can  be  separated  from  that 
due  to  the  target  body.  Consequently  a  cleaned  ISAR  image  of  the  target  body  can  be 
reconstructed  by  using  only  those  bases  associated  with  the  body  motion.  Furthermore, 
robust  Doppler  information  extraction  can  be  achieved  by  applying  the  period  detection 
algorithm  to  the  component  associated  with  rotating  parts  only.  We  have  applied  this 
algorithm  to  both  simulated  data  from  the  radar  prediction  code  Xpatch  and  some  real 
measurement  results  for  helicopters  and  propeller  airplanes  with  good  success.  The 
measured  data  were  provided  by  Dr.  Victor  Chen  of  Naval  Research  Laboratory.  An 
example  demonstrating  the  algorithm  is  shown  in  Fig.  2.1.1.  Fig.  2.1.1(a)  shows  the 


original  ISAR  image  of  a  helicopter.  The  Doppler  smearing  due  to  the  rotating  blades  is 
clearly  visible  and  overlaps  with  the  target  features.  Fig.  2.1.1(b)  shows  the  image  after 
the  adaptive  joint  time-frequency  feature  extraction  procedure  is  applied  to  remove  the 
Doppler  component  due  to  the  rotating  blades.  The  Doppler  smearing  is  significantly 
reduced  and  the  target  features  become  much  more  apparent.  The  removed  Doppler 
smearing  contains  useful  information  on  the  blade  rotation  rate  and  a  period  detection 
algorithm  can  be  applied  to  extract  the  rotation  rate  of  the  blades.  Fig.  2.1.1(c)  shows  the 
autocorrelation  function  versus  dwell  time  of  the  extracted  blade  contribution.  From  this 
plot,  we  can  determine  the  periodicity  of  the  signal  to  be  0.056  s.  The  reciprocal  of  the 
period  is  17.7  rps  and  it  corresponds  to  the  product  of  the  blade  rotation  rate  and  the 
number  of  blades  for  this  target.  This  algorithm  can  also  be  extended  to  deal  with  jet 
engine  modulation  (JEM)  and  is  discussed  next. 

2.2.  Removal  of  Image  Artifacts  due  to  Jet  Engine  Modulation.  Jet  engine 
modulation  (JEM)  is  a  well-known  phenomenon  caused  by  the  high-speed  rotation  of  the 
aircraft  engine.  For  imaging  radars,  the  typical  PRF  is  much  slower  than  the  engine 
rotation  frequency.  Thus  the  resulting  ISAR  image  in  the  frontal  region  of  an  aircraft 
contains  an  aliased  component  along  the  cross  range  dimension,  as  shown  by  the  TIRA 
image  in  Fig.  2.2.1(a).  Such  effect  is  difficult  to  predict  accurately  using  simulation  [6, 
7].  Furthermore,  JEM  lines  are  noise-like  and  can  corrupt  the  geometrical  features  of  the 
target  in  the  ISAR  image.  For  ISAR-based  target  recognition,  it  would  be  useful  to  devise 
an  algorithm  to  separate  the  JEM  lines  from  the  target  image  before  the  subsequent 
classification  process.  In  2.1,  we  demonstrated  the  separation  of  the  rotating  blade 
contribution  from  the  body  image  in  helicopter  data.  However,  JEM  possesses  new 
challenges  due  to  the  high  rotational  rate  of  the  jet  engine  and  additional  electromagnetic 
propagation  effect  through  the  inlet  duct.  We  have  carried  out  JEM  removal  on  TIRA 
data  using  the  AJTF  algorithm.  The  model  we  adopt  assumes  that  the  aircraft  consists  of 
a  slowly  rotating  body  with  a  constant  rotational  velocity  Q.b  and  a  fast  moving  engine 
component  with  a  different  rotational  velocity  Qp.  The  received  radar  return  as  a  function 
of  dwell  time  can  thus  be  written  as: 


(2.2.1) 


^(^d)  =  Li^k  exp[-7— (/?(?o)  +  \  cos(Q^^o)  +  Jk 

k=\  ^ 

+  X  exp[-y  — (/?(?o )  +  cos(Q^r^ )  +  y^  smiQ^^t^ )] 

A=iV>+l  C 

where  N  is  the  total  number  of  point  scatterers  within  one  range  ceU,  of  which  Nb  are  the 
body  scatterers.  Usually  Qp  is  much  greater  than  Qb-  While  the  first  term  can  be 
meaningfully  mapped  into  the  image  plane  of  the  target  via  the  Fomier  transform,  the 
second  term  results  in  serious  Doppler  smearing  across  the  cross  range  domain. 

We  can  also  utilize  the  AJTF  technique  to  separate  the  fast  nooving  part  from  the 
relatively  slow  moving  body.  For  the  component  due  to  target  body  scattering,  the 
Doppler  frequency  is 

/d  +  +  (2.2.2) 

C  c 

while  the  Doppler  frequency  due  to  the  fast  rotating  part  is 
AtiF 

fa  =—^^p[ycos{Qptj^)  +  xsm{Qpt^)]  (2.2.3) 

We  see  that  (2.2.2)  is  a  linear  function  of  dwell  time  while  (2.2.3)  is  a  sinusoidal 
function.  If  we  parameterize  the  signal  by  basis  functions  that  have  linear  Doppler 
frequency  behavior  as  a  function  of  dwell  time,  the  two  signals  can  be  approximately 
separated  by  their  displacement  and  slope  parameters.  We  utilize  the  AJTF  processing 
technique  to  carry  out  the  parameterization.  The  signal  component  due  to  the  target 
scattering  is  reconstructed  by  using  all  the  bases  with  small  displacement  and  small  slope 
parameters.  Fig.  2.2. 1(b)  shows  the  image  after  the  JEM  removal  processing.  Note  that 
the  body  features  are  unveiled  after  the  JEM  removal.  Fig.  2.2.2  shows  the  correlation 
between  the  synthetic  images  and  the  measured  images  from  I'lRA  after  JEM  removal. 
We  observe  that  the  correlation  coefficients  in  the  JEM  region  are  increased  after  we 
remove  the  JEM  interference  from  the  body.  Further  research  is  needed  in  this  area  to 
more  definitively  assess  the  effect  of  JEM  removal  in  the  performance  of  target 
classification  algorithms. 
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Fig.  2.2.1.  TIRA  image  from  a  frontal  view,  (a)  Before  JEM  processing, 
(b)  After  JEM  removal  using  the  AJTF  algorithm. 
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Fig.  2.2.3.  Correlation  between  JTF  and  Xpatch  images 
after  JEM  removal. 

2.3.  Extraction  of  MicroDoppler  Features  from  ISAR  Data  Using  Adaptive  Chirplet 
Representation.  Under  the  assumption  of  rigid  body  motion,  traditional  high-resolution 
ISAR  imaging  is  capable  of  generating  a  (range)-(cross  range)  image  of  the  target.  When 


the  target  contains  non-rigid  body  motion,  the  resulting  radar  image  does  not  preserve  the 
spatial  features  of  the  target.  In  2.1  and  2.2,  we  have  studied  algorithms  to  remove  the 
artifacts  due  to  the  rotating  blades  of  a  helicopter  and  the  jet  engine  modulation  (JEM). 
In  this  work,  we  set  out  to  extract  the  so-called  microDoppler  features  due  to  the  small 
motions  of  a  non-rigid  target  [8]. 

For  this  purpose,  the  chirplet  basis  function  proposed  in  [9]  is  used.  The  radar  signal 
as  a  function  of  dwell  time  t  is  expanded  in  terms  of  N  basis  functions  as  follows: 

s{t)  =  exphOTj  (r  -  tkf  +  {t  - f]  (2.3.1) 

k=\ 

Each  basis  function  is  a  four-parameter  chirplet  with  time  location  tk,  frequency  center /<:, 
time  extent  oCk,  and  chirp  rate  Pk-  The  parameters  for  all  the  individual  chirplets  are  found 
adaptively.  The  chirplet  with  the  highest  energy  is  extracted  first.  We  exhaustively 
search  for  the  maximum  projection  from  the  radar  signal  onto  the  set  of  basis  functions  in 
the  parameter  space.  Once  the  best  basis  is  identified,  we  then  subtract  its  contribution 
from  the  radar  signal  and  continue  to  search  for  the  next  best  basis.  This  process  is 
iterated  until  we  have  a  good  representation  the  original  signal. 

The  radar  data  from  the  Pl-walking  data  collected  using  the  APY-6  radar  has  been 
used  to  test  the  algorithm.  During  the  data  collection,  a  man  walked  toward  the  radar.  At 
least  two  different  motions  were  present  on  the  walking  man:  the  steady  body  motion  and 
the  swinging  arm  motion.  We  first  carry  out  a  coarse  range  alignment  by  correlating  the 
range  profiles.  The  signal  through  a  particular  range  cell  is  then  used  for  JTF  processing. 
Shown  in  Fig.  2.3.1  is  the  spectrogram  of  the  radar  signal  in  a  particular  range  cell 
containing  strong  micro-Doppler.  We  immediately  recognize  the  two  motion  components 
present  in  this  figure.  The  horizontal  trajectory  is  due  to  the  body  motion,  while  the 
sinusoidal  curve  is  due  to  the  arm  swing. 

We  apply  the  JTF  extraction  and  extract  50  chirplets.  After  the  parameterization, 
we  separate  the  body  and  the  arm  returns.  As  can  be  seen  in  Fig.  2.3.1,  the  body  return 
should  consist  of  chirplet  bases  with  both  small  Doppler  frequency  fk  and  small  chirp  rate 
pk-  By  applying  this  criterion  to  the  basis  functions  representing  the  signal,  we  can 
separate  the  body  return  from  the  arm  return.  The  results  are  shown  in  Figs.  2.3.2(a)  and 
2.3.2(b),  respectively.  Notice  that  the  main  features  in  the  two  returns  are  preserved  after 


the  JTF  extraction.  The  Doppler  features  due  to  the  arm  swing  shown  in  Fig.  2.3.2(b)  is 
nearly  periodic.  This  period  can  be  easily  estimated  from  the  arm-only  data  by  taking  the 
autocorrelation  of  the  time  sequence.  The  result  is  shown  in  Fig.  2.3.3(a).  The  peaks  in 
the  autocorrelation  function  indicate  the  period  of  the  motion  and  it  is  estimated  to  be 
0.44  sec.  Note  that  without  using  JTF  analysis,  accurate  detection  of  the  swing  period 
would  be  more  difficult  due  to  the  presence  of  the  body  return.  The  autocorrelation 
calculated  from  the  unfiltered  data  is  shown  in  Fig.  2.3.3(b).  As  we  can  see,  the  peaks  are 
much  less  prominent. 
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Fig.  2.3.1.  Spectrogram  of  the  radar  signal  showing 
micro-Doppler 


Fig.  2.3.2.  Separation  of  the  contribution  from  (a)  the  body  and 
(b)  the  swing  arm  using  the  adaptive  chirplet 
representation. 
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Fig.  2.3.3.  Autocorrelation  to  determine  the  rate  of  arm 
swing,  (a)  after  JTF  extraction  (b)  before  JTF 
extraction. 


3.1.  Three-Dimensional  Motion  Detection  Using  JTF  Algorithm.  One  basic 
assumption  of  existing  motion  compensation  algorithms  is  that  the  target  only  undergoes 
motion  in  a  two-dimensional  plane  (i.e.,  it  has  a  steady  rotational  axis)  during  the  dwell 
interval  needed  to  form  an  image.  From  several  independent  examinations  of  measured 
ISAR  data  sets  recently,  it  was  reported  that  the  presence  of  three-dimensional  motion  is 
quite  detrimental  to  the  formation  of  a  well-focused  image  [10-12].  This  is  also  consistent 
with  our  own  findings  from  the  TIRA  and  MERIC  data.  For  some  image  frames,  we 
found  that  the  motion  compensated  image  and  the  associated  truth  images  are  both  quite 
poor.  An  example  image  is  shown  in  Fig.  3.1.1(a).  When  we  examine  the  motion  data 
from  the  aircraft  instrument,  we  find  that  the  aircraft  motion  is  not  confined  to  a  two- 
dimensional  plane  during  the  imaging  interval.  Fig.  3.1.1(b)  shows  the  plot  of  the 
elevation  versus  azimuth  look  angles  of  the  aircraft  from  the  radar.  This  curve  should  be 
linear  if  the  motion  is  strictly  two-dimensional,  and  in  this  case,  there  is  clearly  three- 
dimensional  motion.  In  general,  we  will  not  be  able  to  form  a  focused  image  using  an 


existing  motion  compensation  assumption,  since  the  assumed  model  is  mismatched  to  the 
actual  target  motion.  The  ultimate  goal  in  this  research  is  to  develop  a  general  motion 
compensation  algorithm  to  handle  targets  with  arbitrary  three-dimensional  motion. 
However,  this  problem  is  quite  challenging,  if  not  impossible  [13].  Instead,  we  address  a 
less  ambitious  problem  of  detecting  the  presence  of  three-dimensional  motion  from  raw 
radar  data.  It  is  hoped  that  the  solution  here  will  be  a  stepping  stone  to  the  ultimate  three- 
dimensional  motion  compensation  problem. 

Allowing  for  arbitrary  three-dimensional  motion  in  space,  we  consider  the  following 
model  as  a  generalization  of  the  two-dimensional  motion  model: 

^  47tf 

E(ta)=  'LA^expl-j - )  +  Zk</>(to  ))]  (3.1.1) 

k=l  C 

where  6  is  the  azimuth  angle  of  the  target  with  respect  to  the  radar,  and  (p  is  the  elevation 
angle.  In  (3.1.1),  it  is  assumed  that  the  translation  motion  has  been  removed  and  that  the 
standard  small-angle,  small  bandwidth  approximations  apply.  This  model  reduces  to  the 
standard  two-dimensional  motion  model  when  6  and  (p  are  linearly  related.  Our 
approaeh  to  the  three-dimensional  motion  detection  problem  is  to  utilize  the  AJTF 
algorithm  to  extract  the  phase  behavior  of  the  radar  data  at  multiple  range  cells.  It  can  be 
shown  that  when  the  target  undergoes  only  two-dimensional  motion  during  the  dwell 
duration,  the  relationship  between  the  phase  extracted  from  one  range  eell  and  that  from 
another  range  eell  should  be  linear.  For  three-dimensional  motion,  the  relationship  is  in 
general  nonlinear.  Therefore,  by  examining  the  linearity  of  the  phase  relationships  from 
different  range  eells,  we  can  distinguish  two-dimensional  motion  from  three-dimensional 
motion.  Fig.  3.1.2  shows  our  results  from  the  TIRA  data.  Fig.  3.1.2(a)  shows  the  degree 
of  three-dimensional  motion  in  the  data  for  20  different  image  frames,  deteeted  by 
applying  our  algorithm  to  the  raw  TIRA  radar  data.  As  a  referenee  for  comparison.  Fig. 
3.1.2(b)  shows  the  degree  of  three-dimensional  motion  for  the  same  20  frames  measured 
using  the  instrumentation  data.  It  ean  be  seen  that  our  algorithm  correctly  detects  where 
significant  three-dimensional  motions  exist.  We  believe  this  detection  algorithm  could  be 
quite  useful  for  determining  the  “good”  imaging  intervals  from  which  focused  images  can 
be  more  readily  generated.  For  targets  that  exhibit  very  chaotic  motions,  such  as  ships  on 


the  ocean,  finding  such  intervals  of  opportunity  may  be  very  critical  for  target 
recognition. 
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Fig.  3.1.1.  Effect  of  non-steady  imaging  plane  on  image  quality,  (a)  A  poorly 
focused  image,  (b)  The  corresponding  target  motion  from  the 
instrumentation  data. 
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(a)  Blind  detection  from  radar  data 
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Fig.  3.1.2.  Blind  detection  of  three-dimensional  motion  from  TIRA  data,  (a) 
Degree  of  three-dimensional  motion  over  20  image  frames  detected 
using  the  proposed  algorithm,  (b)  Degree  of  three-dimensional 
motion  measured  from  the  instrument  data. 


3.2.  3D  ISAR  Image  Reconstruction  of  a  Tai^et  with  Motion  Data.  In  the  previous 
topic,  3D  target  motion  was  treated  as  a  problem  for  ISAR  imaging  in  the  sense  that 
traditional  image  formation  algorithm  with  a  2D  motion  model  cannot  accommodate  the 
situation.  The  algorithm  developed  in  3.1  aUows  such  intervals  to  be  detected  so  that  we 
can  avoid  them  altogether.  In  fact,  the  presence  of  chaotic  target  motion  actually  provides 
an  opportunity  for  3D  ISAR  image  formation,  since  radar  data  is  available  in  a  2D 
angular  aperture.  In  this  topic,  we  try  to  devise  an  algorithm  to  form  a  3D  image  of  a 
target.  We  assume  that  the  motion  data  of  the  target  is  known.  3D  ISAR  imagmg  is 
straightforward  based  on  the  Fourier  transform  of  the  data  in  the  (frequency)-(azimuth)- 
(elevation)  domain.  However,  the  problem  with  the  Fourier  method  is  data  availability.  In 
practice,  the  available  radar  data  is  usually  severely  undersampled  over  the  2D  angular 
aperture.  Therefore,  we  apply  the  adaptive  feature  extraction  (AFE)  algorithm  to  attack 

this  problem 

The  adaptive  feature  extraction  algorithm  is  a  model-based  signal  processing  method 
very  similar  to  the  AJTF  engine  used  for  ISAR  motion  con5)ensation.  The  difference  is 
that  we  use  a  different  class  of  basis  functions,  namely,  the  radar  return  from  a  point 
scatterer  with  unknown  cross  range  positions  but  with  known  motion.  The  two 
orthogonal  cross  range  positions  are  obtained  by  a  search  procedure.  The  strongest  point 
scattorer  is  extracted  when  the  projection  from  the  radar  signal  onto  the  basis  function  is 
maximized.  We  then  subtract  the  strongest  point  scatterer  and  iterate  the  process  for  the 
subsequent,  weaker  point  scatterers. 

We  have  successfully  reconstructed  a  3D  image  of  an  air  target  from  its  high- 
resolution  radar  data  using  this  algorithm  Fig.  3.2.1  shows  the  motion  data  corresponding 
to  multiple  flights  of  the  aircraft.  From  Fig  3.2.1,  we  notice  that  the  motion  data  is  non- 
uniformly  distributed  in  the  5  degree  by  5  degree  (azimuth,  elevation)  window.  After 
applying  the  3D  AFE  algorithm  to  the  available  data,  we  achieve  the  results  as  shown  in 
Fig.  3.2.2.  Fig.  3.2.2(a)  shows  the  extracted  3D  image.  Figs.  3.2.2(b)-(d)  show  the 
projected  top  view,  side  view  and  front  view  of  the  target,  respectively.  In  these  images, 
the  fuselage,  two  wings  and  the  tail  structures  are  correctly  reconstructed.  Through  the 
3D  image  formation,  more  features  are  made  available  than  either  the  2D  ISAR  image  or 
ID  range  profiles.  This  suggests  that  3D  imaging  is  desirable  for  target  identification 


purpose.  This  work  also  provides  an  important  stepping  stone  toward  the  ultimate  goal  of 
blind  3D  image  formation  without  any  motion  data. 


Fig.  3.2.1.  Motion  data 
from  multiple  flights 
of  an  air  target. 
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Fig.  3.2.2.  Imaging  result  with  3D  AFE  applied  to  measurement  data, 
(a)  Extracted  3-D  Image,  (b)  Top  view,  (c)  Side  view,  (d)  Front  view. 


3.3.  Application  of  Genetic  Algorithms  to  ISAR  Imaging.  While  the  AJTF  algorithm 
has  been  utilized  extensively  in  our  work,  a  major  drawback  with  this  method  is  the  high 
computational  complexity  associated  with  the  parameter  search  procedure.  This  problem 
becomes  more  acute  when  the  order  of  the  motion  is  high.  We  have  investigated  the  use 
genetic  algorithms  (GA)  [14]  in  the  search  procedure.  In  contrast  to  conventional 
optimization  methods,  GA  is  a  population-based,  statistical  search  technique.  It  borrows 
such  concepts  as  inheritance  and  mutation  from  the  biological  evolution  process.  As  a 
global  optimization  technique,  GA  is  known  to  be  very  easy  to  implement  and  applicable 
to  many  design  and  inverse  problems  [15].  For  our  problem,  the  objective  function  is 
defined  as  the  projection  magnitude  from  the  radar  signal  onto  the  time-frequency  basis 
functions.  Both  real-coded  and  binary-coded  GA  have  been  implemented  and  the 
performance  has  been  compared  with  results  from  the  exhaustive  search. 

Point  scatterer  simulations  are  first  used  to  test  the  use  of  GA  for  ISAR  motion 
compensation.  Even  for  a  simple  phase  estimation  problem,  the  objective  function  shows 
many  local  maxima,  implying  that  such  problems  would  be  challenging  for  gradient- 
based  local  optimization  techniques.  Fig.  3.3.1(a)  shows  that  the  resulting  phases  from 
both  the  real  and  binary  GA  agree  well  with  the  original  phase  function.  In  Fig.  3.3.1(b), 
we  compare  the  computational  complexity  of  GA  to  exhaustive  search  for  different 
orders  of  motion.  The  exhaustive  search  is  known  to  have  an  exponential  complexity  of 
0(exp(n)).  As  expected,  the  resulting  computation  time  in  logarithmic  scale  shows  up  as 
a  straight  line  in  Fig.  3.3.1(b).  Also  plotted  in  Fig.  3.3.1(b)  are  the  numerical  results  from 
both  binary  and  real-coded  GA.  It  is  observed  that  GA  has  much  lower  complexity  than 
exhaustive  search.  We  also  found  that  real-coded  GA  consistently  achieves  better 
accuracy  than  the  binary-coded  GA.  This  is  because  real-coded  GA  has  the  ability  to 
search  for  any  real  values  within  the  search  range. 

We  have  also  applied  GA  to  the  NATO  TIRA  data.  Fig.  3.3.2(a)  shows  the  ISAR 
image  from  a  particular  time  interval,  along  with  the  spectrogram  of  the  radar  signal  at 
range  cell  64  after  third-order  translation  motion  compensation.  We  see  that  the  point 
scatterers  in  the  image  are  not  focused.  The  trajectory  in  the  spectrogram  is  also  tilted. 
This  means  that  third-order  motion  model  is  not  sufficient  to  describe  the  target  motion. 
After  fourth-order  translation  and  rotational  motion  compensation,  the  result  is  shown  in 
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Fig.  3.3.2(b).  We  see  much  better-focused  point  scatterers  in  the  ISAR  image.  The 
associated  spectrogram  shows  a  straight,  horizontal  trajectory.  While  real  GA  took  45 
seconds  of  computation  time  to  obtain  the  correct  phase,  the  computation  time  for  fourth 
order  exhaustive  search  is  over  50  minutes.  Therefore,  the  time  savings  of  GA  over 
exhaustive  search  becomes  quite  significant  when  the  target  exhibits  highly  irregular 
motion  during  the  imaging  interval. 

4.1.  Clutter  Reduction  for  SAR  Images  Using  Adaptive  Wavelet  Packet  Transform. 

Synthetic  Aperture  Radar  (SAR)  images  of  ground  targets  generally  consist  of  target 
features  and  clutters  from  background  scattering.  In  automatic  target  recognition  (ATR) 
applications,  it  is  desirable  to  remove  the  clutter  from  the  target  image  before  ATR 
processing.  The  standard  way  to  suppress  clutter  is  to  apply  an  appropriate  threshold  level 
to  the  whole  SAR  image.  However,  this  approach  assumes  that  the  target  signal-to-clutter 
ratio  (SCR)  is  large  enough.  Otherwise  this  direct  threshold  approach  results  in  either 
target  feature  loss  or  remnant  clutter  residue.  In  this  work,  we  set  out  to  develop  a 
decluttering  algorithm  to  automatically  extract  the  target  image  from  a  SAR  image  by 
maximizing  the  SCR  using  the  adaptive  wavelet  packet  transform  (AWPT)  [16].  The 
wavelet  packet  basis  is  a  generalization  of  the  conventional  wavelet  basis  [17]  and  has 
been  applied  for  image  compression  [18]  and  moment  matrix  sparsification  [19].  Our 
approach  is  to  transform  the  SAR  image  to  a  new  domain  using  the  wavelet  packet  basis. 
Since  a  typical  target  image  usually  consists  of  point  scatterers  and  more  diffused  region 
features,  the  multi-scaled  wavelet  basis  is  well  suited  to  focus  the  target  image.  Clutter 
image,  on  the  other  hand,  is  statistically  uncorrelated  from  pixel  to  pixel,  and  the 
transformed  clutter  image  under  the  same  set  of  bases  remains  unfocused.  Therefore,  we 
expect  that  the  SCR  can  be  increased  by  transforming  the  original  image  using  an 
appropriately  chosen  set  of  wavelet  packet  basis.  The  cost  function  of  our  AWPT 
algorithm  is  chosen  to  describe  how  well  the  target  signal  is  focused  in  the  transform 
domain.  An  efficient  basis  search  algorithm  is  implemented  to  find  the  best  wavelet 
packet  basis.  Our  algorithm  is  tested  using  the  MSTAR  SAR  data  set  [20].  Fig.  4.1.1(a) 
shows  an  MSTAR  image  in  which  the  target  is  a  ground  vehicle  and  the  clutter  is  due  to 
vegetation.  There  are  several  strong  point  scatters  in  the  front  of  the  vehicle,  but  the 


scattering  from  the  back  part  is  relatively  weak.  Fig.  4.1.1  (b)  shows  the  result  of  applying 
the  direct  thresholding  method  to  the  image.  Fig.  4.1.1(c)  shows  the  decluttered  image  by 
applying  the  AWPT  algorithm.  We  choose  Daubechies  filter  with  order  of  6  as  the 
wavelet  filter.  By  visually  comparing  Figs.  4.1.1(b)  and  (c),  we  note  that  some  crucial 
features  of  the  target  are  kept  in  the  AWPT-processed  image.  In  both  processing 
methods,  there  is  some  target  information  loss.  Fig.  4.1.1(d)  shows  the  signal-to-clutter 
ratio  versus  average  target  image  loss  for  the  two  processing  methods.  It  is  observed  that 
for  a  fixed  target  image  loss  the  AWPT  method  always  achieves  a  higher  SCR  value  than 
the  direct  thresholding  method.  Similar  results  are  obtained  when  the  algorithm  is  applied 
to  other  MSTAR  targets. 


Fig.  4.1.1(a).  SAR  image  of  a  ground 
vehicle  with  clutter. 

Fig.  4.1.1(b).  Clutter  rejection  using 
the  direct  thresholding  method. 
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the  AWPT  approach. 

loss  for  the  two  processing  methods. 

4.2.  Inverse  Scattering  Using  Genetic  Algorithms.  As  an  extension  of  our  research 
into  radar  imaging  and  GA,  we  have  also  carried  out  an  exploratory  study  on  applying 
GA  in  conjunction  with  a  computational  electromagnetics  (GEM)  solver  to  image  targets 
with  strong  multiple  scattering  effects.  Although  this  is  not  directly  related  to  the  JTF 
framework,  radar  imaging  can  be  considered  as  an  approximate,  linearized  version  of  the 
general  inverse  scattering  problem.  Therefore,  our  inverse  scattering  work  provides  us 
with  more  physical  insights  into  the  radar  imaging  problem.  It  is  well  recognized  that 
traditional  imaging  algorithms  suffer  from  resolution  limitation  and  image  artifacts  due  to 
multiple  scattering  phenomena.  Rigorously  solving  the  electromagnetic  inverse  scattering 
problem,  on  the  other  hand,  is  much  more  challenging.  In  this  work,  we  use  GA  together 
with  a  CEM  solver  to  attack  the  two-dimensional  inverse  scattering  problem.  The  method 
of  moements  (MoM)  is  used  for  the  forward  scattering  computation,  while  GA  is  used  as 
a  search  engine  to  minimize  the  difference  between  the  measured  data  and  the  computed 
scattered  fields  from  each  candidate  shape.  We  use  the  binary  bitmap  approach  to 
discretize  the  search  space.  To  constrain  the  problern,  a  geometrical  median  filter  is 
applied  to  create  more  realizable  shapes. 

We  have  applied  our  algorithm  to  actual  measurement  data  from  the  public  Ipswich 
data  set.  Fig.  4.2. 1  (a)  shows  the  shape  and  size  of  three  metallic  Ipswich  objects  selected 
for  inversion,  namely,  the  triangular  cylinder,  the  dihedral  and  the  circular  cavity.  They 
are  labeled  as  Ips009,  IpsOOd  and  IpsOll,  respectively.  We  first  tested  the  inversion 
algorithm  using  MoM-simulated  field  data  as  the  input.  The  search  area  was  chosen  to  be 
15cmxl5cm  for  Ips004  and  12cmxl2cm  for  Ips009  and  IpsOll.  The  number  of  cells 
within  this  area  was  set  to  20x20.  The  reconstructed  results  in  Fig.  4.2.1(b)  show  the  final 
inverted  shapes  of  the  three  objects.  We  observe  that  the  final  shapes  are  in  fairly  good 
agreement  with  the  real  shapes.  The  dihedral  and  the  circular  cavity  contain  strong 
multiple  scattering  and  yet  their  inverted  shapes  closely  resemble  the  correct  objects. 
Results  for  these  targets  were  also  generated  using  the  traditional  imaging  method  and 
they  showed  strong  image  artifacts  due  to  multiple  scattering. 

Next,  we  applied  the  inversion  algorithm  to  the  actual  measured  data.  Fig.  4.2.1(c) 
shows  the  final  reconstruction  shapes.  As  we  can  see,  the  inverted  shape  is  good  for  the 
triangular  cylinder,  which  has  no  multiple  scattering  effects.  For  the  dihedral,  the 


reconstructed  shape  is  not  continuous,  but  is  quite  similar  to  the  real  object.  The  circular 
cavity  shows  the  most  discrepancy  with  the  real  shape.  The  exterior  and  the  opening  of 
the  cavity  are  correctly  inverted,  while  the  interior  part  of  the  cavity  shape  is  not  as 
satisfactory.  Future  work  on  this  topic  should  be  focused  on  finding  ways  to  cut  down  the 
computational  load  such  that  this  method  can  be  extended  to  deal  with  realistic  3D 
targets.  To  accelerate  the  GA  process,  better  geometry  constraints  should  be  applied.  In 
addition,  the  use  of  an  approximate  CEM  solver  such  as  Xpatch  should  be  investigated  to 
handle  large  3D  targets. 


Fig.  4.2.1  Inversion  of  three  Ipswich  objects 

a  Real  shapes  of  Ips009,  Ips004  an  IpsOl  1 
b  Inversion  based  on  MoM-simulated  field 
c  Inversion  based  on  measured  complex  field 


C.  FOLLOW-UP  STATEMENT: 

In  the  past  three  years,  we  have  processed  a  number  of  measurement  data  sets  and 
have  identified  several  important  fundamental  problems  in  radar  imaging.  We  have  also 
developed  a  number  of  JTF-based  algorithms  to  address  these  problems.  We  are  currently 
investigating  the  microDoppler  phenomenon,  which  arises  in  targets  that  violate  the  rigid 
body  model  of  scattering.  MicroDoppler  can  arise  in  situations  where  moving 
components  (e.g.,  rotating  antennas,  spinning  rotor  blades)  exist  on  a  target,  or  when  the 
target  undergoes  strong  flexing  and  vibration  due  to  motion  dynamics.  We  have 
processed  the  “PI  walking”  data  from  the  Navy  APY-6  radar.  We  plan  to  continue  this 
effort  in  a  new  proposed  program  to  apply  JTF  techniques  for  analyzing  and  extracting 
microDoppler  features  and  to  utilize  such  information  to  improve  ATR  performance. 
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Abstract — ^In  this  paper,  three-dimensional  (3D)  inverse  synthetic 
aperture  radar  (ISAR)  image  reconstruction  with  taiown  motion  data 
is  studied.  In  traditional  two-dimensional  (2D)  ISAR  imaging,  a 
2D  point  scatterer  model  is  adequate  to  consider  the  target  rotation 
motion  with  a  fixed  rotational  axis.  However,  target  motions  with  a 
varying  rotational  axis  are  sometimes  encountered  in  real  situations. 
Under  such  cases,  the  use  of  a  3D  point  scatterer  model  is  necessary 
for  3D  ISAR  image  reconstruction.  An  adaptive  feature  detraction 
algorithm  is  proposed  to  reconstruct  the  3D  image  of  a  tmrget  with 
non-uniformly  imdersampled  radar  data  over  the  azimuth  and  elevation 
aperture.  Simulation  results  based  on  actual  motion  data  of  air  targets 
demonstrate  the  effectiveness  of  the  algorithm. 
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1.  INTRODUCTION 

Inverse  synthetic  aperture  radar  (ISAR)  has  been  identified  as  an 
effective  tool  for  target  recognition  [1,  2].  In  the  ISAR  problem,  the 
radar  is  stationary  and  collects  back-scattering  data  firom  a  moving 
target.  The  geometry  of  such  a  problem  is  shown  in  Figure  1,  where 
the  target  rotates  with  azimuth  angle  ip  and  elevation  angle  d  and  the 
x-axis  is  the  radar  line  of  sight  (LOS).  In  traditional  ISAR  imaging, 
we  assume  that  the  target  rotates  with  a  fixed  rotational  axis  during 
the  image  formation  interval.  This  motion  is  termed  two-dimensional 
(2D)  motion,  since  the  target  motion  is  confined  to  a  2D  plane.  Under 
this  case,  only  2D  ISAR  images  can  be  obtained. 


Figure  1.  Geometry  of  an  ISAR  imaging  problem. 

When  the  target  has  a  varying  rotational  axis  diuring  the  imaging 
interval,  the  target  motion  is  termed  three-dimensional  (3D)  motion. 
This  can  occur  on  a  maneuvering  target.  Three  examples  of  3D  target 
motion  are  shown  in  Figure  2.  Figure  2(a)  shows  a  slight  deviation 
from  a  2D  motion.  This  can  occur  when  an  aircraft  undergoes  a 
well-controlled  maneuver  during  the  imaging  interval  [3].  Figure  2(b) 
shows  a  severe  wave-like  3D  motion.  This  kind  of  motion  is  t3T)ical  of 
ship  targets  due  to  ocean  wave  modulation  [4].  Figure  2(c)  shows  the 
target  motion  data  from  multiple  flight  paths.  Although  each  flight 
path  obej^  2D  motion,  the  cmmdative  data  over  the  angular  aperture 
appear  to  be  3D. 


3D  ISAR  image  reconstruction  of  a  target 
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Figure  2.  Possible  3D  motions  encountered  in  real  ISAR  data 
collection,  (a)  Slight  3D  motion  deviating  from  2D  motion,  (b)  Severe 
wave-like  3D  motion,  (c)  Multiple  flight  paths  of  2D  motion. 


For  the  standard  2D  ISAR  imaging  problem,  3D  motion  of  a 
target  is  considered  undesirable  and  its  detrimental  effect  on  the  2D 
IS.^  imaging  result  has  been  analj^d  in  [5].  To  date,  there  is  no 
good  method  to  deal  with  how  to  image  in  the  presence  of  unknown 
3D  motion  [6].  However,  when  the  motion  data  is  known,  there  is 
an  opportimity  to  produce  a  3D  image.  With  the  known  motion 
data,  the  collected  r^ar  data  is  available  in  the  (frequency)- (azimuth)- 
(elevation)  domain.  Based  on  the  3D  point  scatterer  model,  3D  ISAR 
imaging  is  straightforward  based  on  the  Fourier  transform. 

The  problem  for  3D  ISAR  image  formation  with  the  Fourier 
transform  is  data  availability.  Governed  by  the  Nyquist  sampling 
theorem,  the  Fourier  transform  method  requires  data  that  is  dense 
enough  in  the  angular  aperture.  By  examining  Figure  2,  it  is  clem- 
that  the  available  radar  data  in  practice  are  severely  undersampled  over 
the  2D  angular  aperture.  This  makes  the  Fourier  transform  method 
unsuitable  for  the  3D  ISAR  imaging  problem. 

In  this  paper,  we  set  out  to  reconstruct  the  3D  ISAR  image  of  a 
target  assuming  the  3D  motion  of  the  target  is  known.  Similar  topics 
have  been  reported  in  [7-9].  In  [7],  a  mosaic  3D  image  is  produced  from 
many  2D  images  formed  using  a  super-resolution  algorithm.  In  [8],  the 
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target  height  is  reconstructed  using  multiple  cuts  of  elevation  data. 
In  [9],  a  relaxation-based  algorithm  is  used  to  address  the  so-called 
curvilinear  SAR  problem. 

Our  approach  is  to  use  an  adaptive  feature  extraction  algorithm 
based  on  a  general  3D  point  scatterer  model.  It  is  201  extension  of  our 
previous  work  in  [10],  where  only  2D  motion  is  considered.  Here,  we 
consider  both  the  case  of  severe  and  slight  3D  motions.  When  there  is 
severe  3D  motion  of  the  target,  our  goal  is  to  form  a  3D  ISAR  image. 
When  there  is  only  sHght  3D  motion  that  deviates  from  2D  motion, 
om:  goal  is  to  achieve  a  better  2D  ISAR  image. 

This  paper  is  orgemized  as  follows.  In  Section  2,  the  3D  point 
scatterer  model  is  introduced  to  account  for  3D  target  motion.  The 
adaptive  feature  extraction  algorithm  is  then  described  in  Section  3.  In 
Section  4,  we  show  the  advantages  of  this  algorithm  with  some  simple 
point  scatterer  simrdations  and  discuss  the  limitations  of  the  algorithm. 
We  then  apply  the  algorithm  to  reconstruct  3D  images  of  an  aircraft 
from  realistic  ISAR  motion  data.  Conclusions  are  given  in  the  final 
section. 

2.  3D  POINT  SCATTERER  MODEL 

A  2D  point  scatterer  model  is  usually  used  in  traditional  ISAR  imaging 
[1,  2].  In  this  model,  the  target  consists  of  ideal  point  scatterers.  After 
range  compression,  the  radar  data  can  be  expressed  as 

-E'(io)  =  53  +  yMto)] }  (1) 

i=i  '•  ®  ■’ 

where  fo  is  the  radar  center  frequency  and  to  is  the  dwell  time,  x  and 
y  represent  the  target  range  and  cross-range  positions,  respectively. 
The  target  is  assumed  to  consist  of  N3  point  scatterers,  with  the  i*** 
point  scatterer  depicted  by  position  (x,,  yi)  and  strength  crj. 

The  above  model  is  valid  only  if  the  target  rotational  motion  is 
confined  to  a  2D  plane  and  can  thus  be  described  in  terms  of  only  one 
angular  parameter  (p.  When  there  is  3D  motion  of  the  target,  a  more 
general  3D  model  is  required: 

=  53  +  vM^d)  +  Zi&{tD)]\  (2) 

In  the  above  expression,  a  third  coordinate  z  of  the  target  is  included 
to  represent  the  third  dimension  of  the  target  and  another  independent 
angular  motion  paxcuneter  0  is  introduced  to  describe  the  3D  motion. 
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If  we  know  the  target  angular  motion  parameters  (0,  as  a 
function  of  dwell  time,  then  the  radar  data  is  actually  in  the  (azimuth)- 
(elevation)  domain.  Therefore,  we  can  rewrite  equation  (2)  as 

Ei(p,  =  exp  +  ViV  +  (3) 

Based  on  (3),  we  can  see  that  if  both  9  and  (p  are  evenly  sampled  over 
the  2D  angular  aperture,  we  can  use  a  2D  fast  Fourier  transform  with 
respect  to  9  and  (p  to  form  a  3D  radar  image  in  the  range  x  and  cross 
ranges  y  and  z  domain.  If  the  data  is  unevenly  sampled  yet  dense 
enough,  we  can  still  use  the  Fourier  transform  to  process  the  radar 
data.  However,  when  the  available  data  is  highly  undersampled  as 
those  shown  in  Figure  2,  the  Fourier  transform  method  will  not  produce 
any  meaningful  results.  To  overcome  the  tmdersampling  problem,  the 
so-caUed  adaptive  feature  extraction  algorithm  is  proposed  in  the  next 
section. 

3.  ADAPTIVE  FEATURE  EXTRACTION  ALGORITHM 

Adaptive  feature  extraction  (AFE)  is  a  model-based  signal  processing 
method  that  has  been  used  in  a  previous  work  for  2D  ISAR  imaging 
10].  It  is  similar  to  CLEAN  [11]  and  the  matching  pursuit  algorithm 
12].  After  range  alignment,  we  assume  that  the  range  position  is 
resolved  via  range  compression.  The  two  cross  reinge  dimensions  are 
obtained  by  applying  AFE  based  on  the  3D  point  scatterer  model  to  the 
radar  data  within  a  range  bin.  The  basic  idea  is  to  extract  the  strongest 
point  scatterer  first.  Then  the  response  firom  this  point  scatterer  is 
subtracted  firom  the  total  signal.  The  process  is  then  iterated  for  the 
remaining  point  scatterers. 

Within  a  particular  range  bin,  for  every  possible  point  scatterer 
position  (y,  z),  we  construct  a  basis  function  as  the  radar  signal  from  a 
unit  point  scatterer  located  at  position  (y,  z)  using  the  model  depicted 
in  (3) 

Kto)  —  (4) 

To  find  the  strongest  point  scatterer  within  the  range  cell,  we  search 
for  the  basis  function  that  gives  rise  to  the  maximum  projection  from 
the  radar  data  E{tD)  onto  the  basis.  That  is,  we  can  find  the  position 
of  the  strongest  point  scatterer  (y^,  Zm)  as 

{ym,  Zrn}  =  arg  max  ]<  hito)  >\ 


(5) 
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where  the  projection  is  defined  as 

<  E{tD),h(tD)  >=  J  E{tD)h*{tD)dtD  (6) 

After  we  finri  the  position  of  the  strongest  point  scatterer,  the  strength 
am  of  that  point  scatterer  is  simply  the  projection  from  the  radar  data 
onto  the  chosen  basis  function: 

Cm  =<  E{tD),  >  (7) 

Once  we  have  found  the  contribution  from  the  strongest  point  scatterer 
to  the  radar  signal,  it  is  denoted  as 

Em  =  (8) 

We  then  subtract  Em  from  E  to  find  the  residual  signal: 

Em+l  —  E  Em  (9) 

Finally,  this  process  is  iterated  to  find  the  subsequent  point  scatterers 
one  at  a  time  until  the  energy  of  the  residual  signal  falls  below  a  preset 
threshold. 

To  summarize,  the  steps  in  the  adaptive  feature  extraction 
algorithm  are  as  follows: 

Step  1.  Set  up  the  searching  space  for  both  y  and  z. 

Step  2.  Search  for  metximum  projection  from  the  present  signal 
onto  the  searching  space.  The  two  cross  range  positions  and  the 
strength  of  the  strongest  point  scatterer  within  the  range  bin  are 
determined. 

Step  3i  Remove  the  response  of  the  point  scatterer  from  the  radar 
signal. 

Step  4.  Repeat  steps  2  and  3  until  the  energy  of  the  residual  signal 
falls  below  a  predefined  threshold. 

Step  5.  Repeat  the  above  procedures  for  other  range  bins. 

For  this  work,  the  exhaustive  search  method  is  used  in  step  2. 
It  is  guaranteed  to  produce  a  global  maximum  for  the  projection.  To 
further  save  computation  time,  we  have  also  explored  the  use  of  genetic 
algorithm  instead  of  the  exhaustive  search  for  more  eflSicient  global 
optimization  [13]. 

4.  RESULTS 

To  demonstrate  the  use  of  AFE  for  ISAR  image  formation  of  a  target 
with  3D  motion,  we  first  compare  the  result  against  that  from  the 
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Figure  3.  Original  9  point  scatterers  within  in  one  range  cell,  (a) 
(y,  z)  distribution,  (b)  y  projection,  (c)  z  projection. 


Fourier  transform.  We  then  study  the  resolution  and  noise  sensitivity 
of  the  algorithm.  After  establishing  the  advantages  and  limitations  of 
the  AFE  algorithm,  we  reconstruct  3D  images  from  some  simiilated 
radar  data  based  on  real  motion  data  from  an  air  target. 

4.1.  Testing  on  the  AFE  Algorithm 

To  test  the  algorithm,  we  consider  radar  data  from  a  fixed  range  cell 
and  study  how  to  resolve  the  two  cross  range  dimensions  using  AFE.  A 
target  consisting  of  nine  point  scatterers  of  equal  strength  in  the  (y,  z) 
plane  at  a  fixed  range  is  shown  in  Figure  3(a).  The  two  projections  of 
the  point  scatterers  along  the  y  and  the  z  axes  are  shown  in  Figures 
3(b)  and  3(c)  respectively.  With  the  point  scatterer  model  and  assumed 
3D  motion,  we  generate  simulated  radar  data  based  on  equation  (2). 
The  center  frequency  is  10  GHz  and  the  bandwidth  is  1  GHz  in  our 
simulation.  We  then  apply  the  AFE  algorithm  to  the  simulated  radar 
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Figure  4.  Reconstruction  result  using  AFE.  (a)  Motion  data  in  a 
3  deg  by  3  deg  aperture,  (b)  (y,  z)  image,  (c)  y  projection,  (d)  z 
projection. 


data. 

The  first  case  considered  is  severe  3D  motion  as  shown  by  the 
azimuth-elevation  data  in  Figure  4(a) .  In  this  figure,  256  data  points 
over  a  3-degree  azimuth  by  3-degree  elevation  aperture  are  randomly 
distributed.  The  AFE  reconstructed  image  in  the  (y,  z)  plane  is 
illustrated  in  Figure  4(b).  The  two  projections  along  the  y  and  the 
z  axes  are  shown  in  Figures  4(c)  and  4(d)  respectively.  We  can  see 
that  they  agree  well  with  the  reference  images  in  Figures  3(a)'-3(c). 
The  small  errors  are  due  to  finite  position  resolution  in  the  searching 
space. 

•Next,  we  generate  the  image  via  the  Fourier  transform.  The  result 
of  the  (y,  z)  image  is  shown  in  Figure  5.  It  is  obtained  by  a  brute- 
force  Fourier  transform  over  the  non-uniform  angular  aperture.  The 
sidelobes  are  so  high  that  the  positions  of  the  point  scatterers  cannot 
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Figure  5.  Reconstruction  result  using  the  Fourier  transform,  (a) 
(y,  z)  image,  (b)  y  projection,  (c)  z  projection. 


be  easily  recognized  in  the  (y,  z)  image  or  the  two  projections.  This 
is  due  to  the  kurge  aliasing  effect  associated  with  the  undersampled 
data.  The  Nyquist  sampling  in  angle  is  about  0.07  degree,  while  the 
actual  average  sampling  of  the  data  is  about  0.19  degree.  We  cannot 
e3q)ect  a  good  image  to  be  obtained  by  the  Fourier  transform  method. 
Therefore,  AFE  is  a  much  better  choice  than  the  Fourier  transform  to 
process  the  undersampled  radar  data  at  hand. 

The  second  case  considered  is  slight  3D  motion  as  shown  in  Figure 
6(a).  The  motion  parameters  in  this  figmre  are  taken  from  the  real 
motion  data  of  an  in-flight  aircrafb.  The  azimuth  and  elevation  angles 
are  obtained  through  coordinate  transformation  of  the  original  roll, 
jraw  and  pitch  motions  of  an  instrumented  mrcraft.  As  we  can  see, 
there  is  a  slight  3D  motion  that  deviates  from  the  idealized  2D  motion, 
which  should  be  a  line  in  the  9-(p  aperture.  In  this  case,  we  apply  AFE 
based  on  both  the  correct  3D  motion  model  in  (2)  and  the  approximate 
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Figure  6.  AFE  reconstruction  from  data  with  a  slight  3D  motion,  (a) 
Motion  data,  (b)  {y,  z)  image,  (c)  y  projection,  (d)  z  projection. 


2D  motion  model  in  (1).  In  the  former  case,  the  AFE  search  is  over  the 
(y,  z)  space,  while  in  the  latter  case,  the  search  is  a  one-dimensional 
one  over  the  y  dimension  only.  Figure  6(b)  shows  the  (y,  z)  image 
resulting  from  the  3D  motion  model.  The  two  projections  along  the  y 
and  the  z  axes  are  shown  in  Figures  6(c)  and  6(d).  At  first  glance,  the 
image  in  Figure  6(b)  does  not  look  like  the  original  image  in  Figure  3(a) 
at  all.  However,  we  can  see  that  the  y  projection  in  Figure  6(c)  agrees 
fairly  well  with  the  original  y  projection  in  Figure  3(b).  On  the  other 
hand,  the  z  projection  in  Figure  6(d)  is  very  different  from  the  original 
z  projection  in  Figure  3(c).  The  reason  for  the  poor  performance  is 
that  the  variation  in  is  only  one-fifth  of  the  variation  in  6,  which 
makes  the  image  resolution  in  z  much  lower  than  the  image  resolution 
in  y.  The  resolution  of  the  AFE  algorithm  will  be  discussed  in  more 
detail  in  4.2. 

Figure  7  shows  the  AFE  result  based  on  the  2D  motion  model. 
Only  the  y  projection  is  available  from  the  one-dimensional  AFE 
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Figure  7.  AFE  reconstruction  from  data  with  a  slight  3D  motion 
based  on  a  2D  motion  model. 

search.  The  y  positions  from  this  figure  do  not  refiect  the  true  y 
positions  of  the  9  point  scatterers  shown  in  Figure  3(b).  The  mistake 
here  is  due  to  model  mismatch.  Although  the  variation  in  the  elevation 
angle  $  is  only  one-fifth  of  the  variation  in  the  azimuth  angle  (p,  this 
slight  3D  motion  cannot  be  ignored  in  processing  the  radar  data. 
Therefore,  the  full  3D  motion  model  in  conjunction  with  the  AFE 
algorithm  is  needed  for  image  formation  in  the  presence  of  3D  motio33s. 

4.2.  Resolution  and  Sensitivity  of  the  AFE  Algorithm 

As  we  saw  in  the  last  section,  even  though  the  AFE  algorithm  can 
overcome  the  undersampling  problem,  the  resulting  image  resolution 
is  still  controlled  by  the  apertmre  size.  To  further  illustrate  this,  the  3 
degree  by  3  degree  apertmre  in  Figure  4(a)  is  scaled  to  a  0.3  degree  by  3 
degree  aperture  in  Figure  8(a),  while  stiU  populating  the  aperture  with 
256  data  samples.  The  same  9  point  scatterers  and  radar  parameters 
are  used  to  simulate  the  radar  data.  Only  the  motion  data  is  different. 
The  AFE  reconstruction  results  are  shown  in  Figmes  8(b)-8(d).  As 
we  can  see  by  comparing  Figures  8  to  3,  the  image  resolution  in  the 
y  direction  remains  essentially  unchanged,  while  the  image  resolution 
in  the  z  dimension  gets  much  worse.  Therefore,  the  AFE  resolution 
in  position  is  inversely  proportional  to  the  angular  aperture.  This  is 
exactly  the  same  as  the  Fourier  transform. 

Like  many  other  model-based  signal  processing  techniques,  the 
adaptive  feature  extraction  algorithm  is  sensitive  to  data  error.  We 
study  the  algorithm  sensitivity  by  adding  some  noise  to  the  data. 
Two  t3q>es  of  errors  are  considered.  The  first  is  the  radar  data  error. 
Figure  9(a)  shows  the  reconstructed  image  when  the  radar  data  is 


Figure  8.  Reconstruction  result  using  AFE.  (a)  Motion  data  in  a  3 
deg  by  0.3  deg  aperture,  (b)  (y,  z)  image,  (c)  y  projection,  (d)  z 
projection. 


Figure  9.  Reconstruction  using  noisy  data,  (a)  — 20dB  white  noise 
in  the  radar  data,  (b)  1%  randomness  in  the  angular  positions. 
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Figure  10.  Point  scatterer  representation  of  an  aircraft,  (a)  3D  point 
scatterer  representation,  (b)  Top  view,  (c)  Side  view,  (d)  fVont  view. 


contaminated  with  —20  dB  of  white  noise.  The  second  source  of  error 
is  the  motion  parameter  error.  Figure  9(b)  shows  that  the  result  with 
1%  random  error  in  the  motion  pareuneters.  Our  experience  shows 
that  the  algorithm  is  more  sensitive  to  error  in  the  motion  parameters 
than  that  in  the  radar  data.  To  alleviate  the  angular  motion  error,  we 
suggest  preprocessing  the  motion  data  according  to  a  priori  knowledge 
of  the  motion.  For  example,  we  can  limit  the  motion  data  to  low- 
order  polynomial  functions  of  dwell  time  if  we  believe  the  motion  of 
the  target  is  smooth. 

4.3.  Aircraft  Imaging  with  Real  3D  Motion  Data 

A  full  simulation  is  done  with  point  scatterers  from  an  aircraft  model. 
401-point  scatterers  with  equal  strength  are  used  to  model  oin:  example 
aircraft  as  shown  iu  Figure  10(a).  The  three  projections  of  the  original 
target  corresponding  to  the  top  view,  side  view,  and  front  view  of 
the  aircraft  are  shown  in  Figures  10(b)-10(d).  The  3D  motion  data 


Figure  12.  3D  image  reconstructed  by  AFE  based  on  the  3D  motion 
model,  (a)  3D  view  of  the  top  600  extracted  point  scatterers  (b)  Top 
view,  (c)  Side  view,  (d)  Front  view. 
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in  Figure  11  shows  a  5  degree  by  5  degree  angular  window  in  both 
azimuth  and  elevation.  They  are  from  the  actual  flight  paths  of  an 
aircraft.  To  simulate  the  radar  data,  the  same  radar  parameters  in 
the  previous  examples  are  used.  Our  goal  is  to  reconstruct  a  3D 
image  from  the  range  profiles  simulated  at  these  angles.  Again,  the 
highly  non-uniform,  undersampled  distribution  of  the  radar  data  in 
the  angular  aperture  requires  the  use  of  the  3D  motion  model  and  the 
AFE  algorithm  in  order  to  reconstruct  a  3D  image. 

For  3D  image  construction,  we  need  to  find  the  positions  and 
strengths  of  the  point  scatterers  in  3D  space.  The  range  location  Xs 
is  directly  obtained  from  the  range  profiles.  Then  for  every  range  cell 
within  the  range  profile,  a  collection  of  point  scatterers  with  two  cross- 
range  locations  (ys,  Zs)  and  point  scatterer  strength  Cs  axe  extracted 
from  the  AFE  iterations.  The  computation  time  is  about  15  minutes 
on  a  PC  with  a  P  III  750  MHz  processor.  The  result  is  an  image 
o’s(a^s)  Vs,  Zs)i  representing  the  scattering  strength  of  scatterers  in  3D 
space.  We  sort  the  scattarers  according  to  their  strengths  and  keep  the 
top  600  point  scatterers.  The  positions  of  those  point  scatterers  are 
shown  in  Figure  12(a).  The  top  view,  side  view,  and  front  view  of  the 
reconstructed  3D  image  are  shown  in  Figures  12(b)-12(d),  respectively. 
They  agree  very  well  with  the  three  projections  from  the  original  3D 
point  scatterers  in  Figure  10.  Such  features  as  the  fuselage,  the  wings, 
and  the  vertical  tail  fin  are  correctly  constructed.  Our  algorithm  has 
also  been  tested  on  real  radar  data  of  an  eiir  target. 

5.  CONCLUSIONS 

The  problem  of  3D  ISAR  imaging  of  a  target  with  known  motion 
data  is  discussed  in  this  paper.  The  3D  motion  model  is  necessary  to 
form  a  3D  image  of  the  target  with  3D  motion.  The  adaptive  feature 
extraction  algorithm  along  with  the  3D  motion  model  is  proposed  to 
process  non-uniform,  undersampl^  reidar  data  collected  over  a  2D 
angular  apertme.  The  advantages  of  the  AFE  algorithm  are  clearly 
demonstrated  when  compared  to  either  the  Fourier  transform  or  the 
AFE  based  on  2D  motion  model.  We  find  that  the  resolution  of  the 
AFE  image  is  limited  by  the  aperture  size  in  the  corresponding  motion 
direction.  The  sensitivity  of  the  method  is  also  discussed.  Simulation 
with  a  3D  point  scatterer  model  of  an  aircraft  shows  that  the  resulting 
3D  images  carry  much  more  information  than  the  range  profiOies  or  the 
2D  ISAR  images  alone.  One  advantage  of  the  algorithm  is  that  the 
data  collection  scheme  is  very  flexible.  Therefore,  it  is  a  viable  tool  for 
target  feature  extraction  when  the  3D  motion  data  associated  with  the 
target  are  available. 
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The  significance  of  this  work  is  twofold.  First,  when  there  is 
known  target  motion,  3D  motion  is  taken  into  accormt  based  on  the 
3D  motion  model.  It  is  not  treated  as  an  error  term  as  in  traditional 
2D  ISAR  imaging.  Second,  when  there  is  unknown  target  motion,  this 
work  provides  us  with  an  intermediate  step  toward  blind  3D  ISAR 
image  formation.  Instead  of  doing  the  more  diallenging  3D  motion 
compensation  directly,  we  might  first  estimate  the  3D  motion  fi'om  the 
radar  data.  We  can  then  use  this  method  to  form  a  3D  image. 
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Abstract.  We  present  an  algorithm  to  detect  the  presence  of  3D  target  motion  from 
ISAR  data.  Based  on  the  3D  point  scatterer  model,  we  first  examine  the  effect  of  3D 
motion  on  ISAR  imaging.  It  is  shown  that  existing  motion  compensation  algorithms 
cannot  properly  focus  targets  exhibiting  3D  motion  during  the  imaging  interval.  An 
algorithm  is  then  derived  to  blindly  detect  the  degree  of  3D  target  motion  from  raw  radar 
data.  It  is  based  on  measuring  the  linearity  of  phases  between  two  or  more  point  scatterers 
on  the  target.  The  phase  estimation  is  implemented  using  the  adaptive  joint  time- 
frequency  technique.  Examples  are  provided  to  demonstrate  the  effectiveness  of  the  3D 
motion  detection  algorithm  with  both  simulation  and  real  ISAR  data.  The  detection 
results  are  corroborated  with  the  truth  motion  data  from  on-board  motion  sensors  and 
correlated  with  the  resulting  ISAR  images. 
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1.  Introduction 


High-resolution  inverse  synthetic  aperture  radar  (ISAR)  imaging  is  regarded  as  an 
effective  tool  in  automatic  target  recognition  [1-2].  Ideally,  the  desired  target  motion  is 
uniform  rotation  without  translational  motion,  under  which  a  two-dimensional  (2D) 
Fourier  transform  brings  the  radar  data  in  the  (frequency)-(dwell  time)  domain  into  the 
(range)-(Doppler  frequency)  domain.  Otherwise,  motion  compensation  is  needed  as  an 
intermediate  step  to  form  a  focused  ISAR  image. 

Since  target  motion  can  always  be  decomposed  into  translational  motion  and 
rotational  motion,  a  typical  motion  compensation  algorithm  consists  of  two  steps.  First,  a 
point  on  the  target  is  focused  through  translational  motion  compensation.  When  there  is 
non-uniform  rotational  motion,  other  points  on  the  target  are  not  necessarily  focused. 
Rotational  motion  compensation  is  then  applied  to  focus  these  other  points.  Existing 
motion  compensation  algorithms  usually  assume  that  the  rotational  motion  of  a  target  is 
confined  to  a  2D  plane  during  the  dwell  duration  [1-7],  We  shall  use  the  term  2D  motion 
to  refer  to  target  rotational  motion  of  this  type.  Under  the  2D  motion  assumption, 
rotational  compensation  of  a  second  point  on  the  target  will  focus  the  whole  target. 
When  there  is  3D  motion,  i.e.,  when  the  rotational  motion  is  not  confined  to  a  2D  plane, 
rotational  compensation  of  a  second  point  cannot  focus  the  whole  target. 

Recently,  several  independent  research  groups  have  reported  that,  for  aircraft 
undergoing  fast  maneuvers  or  ships  on  rough  seas,  the  motion  of  a  target  may  be  highly 
chaotic  and  does  not  always  obey  the  2D  motion  model  [8-10].  As  a  result,  the  image 
formed  using  the  standard  motion  compensation  algorithms  is  blurred.  In  [8]  and  [9],  the 
effect  of  3D  motion  on  ISAR  imaging  is  discussed.  However,  target  motions  are  assumed 
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to  be  known  from  other  auxiliary  sensor  data  that  are  usually  not  accessible  in  real 
operational  environment.  In  [10],  the  imaging  interval  is  adaptively  chosen  based  on  the 
resolved  target  feature  in  the  radar  image  to  overcome  the  3D  motion  issue.  It  requires 
sound  knowledge  of  the  target  under  consideration,  which  is  often  not  known  to  the  end 
users  of  ISAR  data. 

The  objective  of  this  paper  is  to  develop  an  algorithm  to  detect  the  presence  of  3D 
motions  during  the  imaging  interval  from  ISAR  data.  Based  on  the  3D  point  scatterer 
model,  we  first  examine  the  effect  of  3D  motion  on  existing  imaging  algorithms.  We  then 
develop  an  algorithm  to  blindly  detect  the  existence  of  3D  motion.  For  this  purpose,  only 
the  estimation  of  phases  of  several  prominent  point  scatterers  is  needed.  It  can  be 
accomplished  by  the  joint  time-frequency  analysis  [6].  With  the  detection  algorithm,  we 
have  the  ability  to  distinguish  the  time  intervals  when  the  target  undergoes  smooth  2D 
motion  from  those  containing  more  chaotic  3D  motion.  As  a  result,  the  good  imaging 
intervals  where  focused  images  are  more  easily  formed  can  be  automatically  determined. 

The  paper  is  organized  as  follows.  First,  the  ISAR  imaging  problem  is  formulated 
in  terms  of  a  point  scatterer  model  in  Section  2.  In  Section  3,  the  2D  motion  assumption 
in  existing  motion  compensation  algorithms  is  analyzed.  We  show  the  reason  why  3D 
motion  is  a  problem  for  ISAR  imaging.  Section  4  discusses  the  3D  motion  detection 
algorithm  in  detail.  Examples  from  both  simulation  and  measurement  data  are  presented 
in  Section  5.  The  conclusions  are  given  in  the  last  section. 

2.  2D  and  3D  Motion  Models 

The  standard  model  used  in  ISAR  processing  is  the  point  scatterer  model  given  as 
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(1) 


E(f,tD)  =  T.^i  ’  yi )  exp  {-j — [K^d  ) + ^;  +  yi<p(tD )] } 

where /is  the  radar  frequency  and  to  is  the  dwell  time.  The  radar  echo  data  E(f,tD)  is  in 
the  (frequency)-(dwell  time)  domain,  x  and  y  represent  the  target  range  and  cross-range 
positions,  respectively.  The  target  consists  of  Ns  point  scatterers,  with  the  i**'  point 
scatterer  depicted  by  position  (x„y,)  and  strength  o(Xi,yi).  The  target  motion  includes  both 
the  translational  motion  described  by  r(tD)  and  the  rotational  motion  described  by  ^to)- 
When  there  is  no  translational  motion  and  the  rotational  motion  is  uniform,  it  is  seen  that 
a  2D  Fourier  transform  brings  the  radar  data  EifJo)  into  a  radar  image  o(x,y).  Otherwise, 
motion  compensation  is  a  critical  step  in  IS  AR  imaging. 

The  above  model  is  what  we  call  a  2D  problem  since  the  target  rotational  motion 
is  confined  to  a  2D  plane  and  describable  in  terms  of  only  one  angular  parameter 
When  there  is  3D  motion  of  the  target,  a  more  general  3D  model  is  required; 

Eif,tD)  =  ^^i(Xi^yi^  Zi )  exp  {-  7  —  )  +  x,  +  y,  ??( )] }  (2) 

In  the  above  expression,  a  third  coordinate  z  of  the  target  is  included  to  represent  the  3D 
target  and  another  independent  angular  motion  parameter  $  is  introduced  to  describe  the 
3D  rotational  motion  (see  Fig.  1).  It  is  possible  to  perform  3D  target  imaging  if  the  target 
motion  is  known  exactly  [11],  [12].  In  practical  ISAR  scenarios,  however,  we  have  no 
access  to  the  target  motion.  Our  objectives  here  are  to  examine  the  effect  of  3D  motion 
on  ISAR  imaging  and  devise  an  algorithm  to  detect  the  presence  of  3D  motion  from  the 
radar  data  itself,  i.e.,  without  any  additional  knowledge  of  the  target  motion. 


4 


3.  Problem  of  Existing  Motion  Compensation  Algorithms  with  3D  Target  Motion 

First,  we  show  that  the  more  general  3D  model  degenerates  into  the  2D  model 
under  two  conditions.  The  first  case  is  when  there  is  a  linear  relationship  between  (p  and 
9,  i.e., 

d{t^)  =  b(p(,t^)  (3) 

This  allows  us  to  cast  equation  (2)  into  the  form 

Ns  A^TTF 

=  (^/  >  yi  ’  Zi )  exp  {-;■ — [r(fo )  +  ^,  +  ( )] }  (4) 

Comparing  (4)  with  (1),  we  see  that  if  we  define  j,-  =  >»,•  +bzi ,  then  the  rotational  motion 

t 

is  in  fact  a  two-dimensional  one  and  the  resulting  2D  image  (T(  xi ,  yi  )  is  the  projection 
from  the  3D  target  trf  xi,yi,Zi )  onto  the  2D  motion  plane. 

The  second  case  is  when  the  z  dimension  of  the  target  is  so  small  that  the  third 
phase  term  in  (2)  can  be  neglected.  For  example,  suppose  a  radar  operates  at  a  frequency 
of  10  GHz  and  the  9  variation  is  limited  to  0.5  degree.  If  the  target  thickness  in  the  z- 
dimension  is  less  than  0.2  m,  then  the  third  phase  term  is  less  than  7c/4  and  the  2D  model 
is  adequate. 

From  the  above  consideration,  we  see  that  the  2D  model  is  applicable  if  either  the 
rotational  motion  is  2D  or  the  target  is  of  2D  in  extent.  When  there  exists  3D  motion  on  a 
full  3D  target,  any  motion  compensation  algorithms  based  on  the  2D  model  is  not 
expected  to  focus  the  target  well.  We  will  now  examine  this  issue  in  more  detail.  Since 
the  translational  motion  compensation  is  independent  of  the  models  in  (1)  and  (2),  only 
the  rotational  motion  compensation  needs  to  be  investigated.  With  2D  rotational  motion 
present,  the  phase  of  a  point  scatterer  i  due  to  the  rotational  motion  is 
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PiitD)  =  yi<p(tD) 


(5) 


Here,  the  constant  4^/c  has  been  suppressed  for  simplicity.  As  we  can  see  from  (5),  the 
phases  of  all  the  point  scatterers  are  linearly  related  to  each  other  (through  the  ratio  of 
their  cross  range  positions).  If  we  make  one  of  the  phases  a  linear  function  of  time,  then 
all  the  phases  are  linearized  simultaneously,  and  the  whole  target  can  be  focused  after  the 
Fourier  transform.  This  is  the  basis  of  most  2D  rotational  motion  compensation 
algorithms  based  on  the  point  scatterer  model  [3-6].  This  concept  is  illustrated  in  Figure 
3.  Figure  3(a)  shows  the  phase  functions  of  two  point  scatterers  under  2D  rotational 
motion.  Figure  3(b)  shows  that  both  points  can  be  made  linear  functions  of  time  after  we 
force  one  of  them  to  be  a  linear  function. 

With  3D  motion,  the  phase  of  a  point  scatterer  due  to  the  rotational  motion  is 

Pi(tD)  =  yi9(^D)  +  ZiO(tD)  (6) 

In  this  case,  the  phases  of  the  point  scatterers  are  no  longer  linearly  related.  If  we  make 
one  of  the  phases  a  linear  function  of  time,  the  phases  of  the  other  point  scatterers  are  not 
automatically  made  linear  functions  of  time,  as  was  the  case  of  2D  motion.  Figure  3(c) 
shows  the  phase  functions  of  two  point  scatterers  with  3D  motion.  As  we  can  see  from 
Figure  3(d),  after  one  point  is  forced  to  be  of  linear  phase,  the  phase  of  the  other  point 
remains  nonlinear. 

Figure  4  illustrates  some  simulation  results  of  the  effect  of  the  rotational  motion 
compensation  based  on  the  model  in  (1)  on  the  final  images  under  2D  and  3D  target 
motion.  The  adaptive  joint  time-frequency  (AJTF)  algorithm  reported  in  [6]  is  used  for 
motion  compensation.  Ten  points  in  3D  space  are  used  to  simulate  the  radar  data.  Figure 
4(a)  shows  an  assumed  2D  rotational  motion.  Figure  4(b)  shows  the  image  after  the 
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translational  motion  compensation.  The  image  shows  one  point  being  focused  in  range 
cell  25  while  other  points  are  unfocused  due  to  the  rotational  motion.  Figure  4(c)  shows 
the  image  after  the  2D  rotational  motion  compensation  in  which  a  point  scatterer  in  range 
cell  57  is  selected  for  focusing.  All  the  point  scatterers  are  focused  in  the  image.  The 
situation  with  an  assumed  3D  target  motion  is  shown  in  Figures  4(d)-4(f).  Figure  4(d) 
shows  the  assumed  3D  motion.  Figure  4(e)  shows  the  image  after  translational  motion 
compensation.  Figure  4(f)  shows  the  final  image  after  the  2D  rotational  motion 
compensation.  The  two  points  in  range  cells  25  and  57  are  focused,  as  expected.  Another 
point  scatterer  in  range  cell  99  is  also  focused  as  it  happens  to  be  in  the  same  2D  motion 
plane  as  the  point  scatterer  in  range  cell  57.  As  we  can  see,  it  is  not  possible  to  foeus  all 
the  points  simultaneously  with  an  existing  algorithm  based  on  the  2D  motion  model. 

4.  3D  Motion  Detection  Algorithm 

Since  existing  motion  compensation  algorithms  cannot  handle  3D  target  motion, 
it  is  desirable  to  develop  a  general  compensation  algorithm  that  can  accommodate  3D 
motion.  However,  this  is  a  difficult  task  (see  [13]  for  background  on  this  problem)  and 
outside  the  scope  of  this  work.  Our  goal  here  is  to  develop  an  algorithm  to  detect  the 
presence  of  3D  motion  from  radar  data.  If  we  can  reliably  detect  those  time  intervals 
where  2D  target  motions  are  predominant,  we  can  use  the  existing  2D  motion 
compensation  algorithms  to  form  well-focused  ISAR  images. 

As  discussed  in  the  last  section,  2D  motion  can  be  represented  by  a  linear 
relationship  between  6  and  Therefore,  we  set  out  to  detect  the  existence  of  a  nonlinear 
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relationship  between  6  and  (p  in  our  3D  motion  detection  algorithm.  First,  we  write  the 
relationship  between  ^and  ^into  a  linear  and  a  nonlinear  part  as  follows: 

0{to)  =  b(p{tD)  +  m(tj;)  (7) 

where  h  is  the  linear  constant  and  mito)  is  the  nonlinear  part  which  indicates  deviation 
from  2D  target  motion,  or  the  degree  of  3D  motion.  Next  we  try  to  gather  target  motion 
information  by  analyzing  the  phases  of  two  point  scatterers  on  the  target.  Let  us  write  the 
relationship  between  the  phase  functions  Pi  and  Pi  of  two  point  scatterers  as: 

P,(rJ  =  aP,(r^)  +  n(r^)  (8) 

The  relationship  is  again  decomposed  into  the  linear  part,  where  a  is  the  linear  constant, 
and  the  nonlinear  part  nito).  Our  goal  is  to  derive  a  relationship  between  m(tu)  and  nitu) 
so  that  the  presence  of  m  can  be  detected  by  observing  n. 

After  the  standard  translational  motion  compensation,  the  time-varying  phase  of  a 
point  scatter  is  in  the  form  of 

/^.(ro)  =  Ay,^(ro)  +  Az,^(fD)  (9) 

where  Ayi  and  Azi  are  differential  positions  of  point  scatterer  i  relative  to  the  reference 
point  chosen  during  translational  motion  compensation.  Substituting  (7)  into  (9)  and  then 
evaluating  (9)  at  point  scatterers  1  and  2,  we  have 

P,  (r^ )  =  (Ay,  +  bAz^  )(p{ta )  +  Az,m(h )  ( 

Pj  (t^ )  =  (Ayj  -I-  bAZi  )  (10b) 

We  next  substitute  (10)  into  (8),  which  leads  to 

a[(  Ay,  +  bAzi  )(p(tjj )  +  Az^mit^ )]  +  n(ro  )  =  (Ay^+  bAz^  )<p(tjj )  +  Azipiitj) )  (11) 
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Notice  that  if  there  is  only  2D  motion,  then  the  phases  of  the  two  point  scatterers  must  be 
linear.  This  means  if  m=0,  then  n=0.  By  using  this  fact  and  equating  the  coefficients  of 
cpitj))  in  (1 1),  the  constant  a  can  be  derived: 


a  = 


Ay,  +  feAz, 

By  substituting  (12)  into  (1 1),  we  finally  arrive  at 
^  ^  Ay,  +  Z^A^i  .  . 

^  «(^D ) 

Az2Ay,  -  Ay^Az, 


(12) 


(13) 


Equation  (13)  states  that  once  the  nonlinear  phase  term  n  is  known,  it  is  proportional  to 
nonlinear  motion  m.  Therefore,  the  steps  to  determine  the  degree  of  3D  target  motion  are 
as  follows.  First,  we  extract  the  phases  of  two  point  scatterers  from  the  radar  data.  Next 
we  find  the  nonlinear  phase  function  n  using  a  minimum  least  squares  fit  of  equation  (8). 
Once  n  is  known,  we  use  equation  (13)  to  decide  on  the  degree  of  3D  motion.  The 
remaining  issues  are:  (i)  how  to  determine  the  phase  functions  of  the  point  scatterers,  (ii) 
how  to  define  the  degree  of  nonlinearity  and  the  degree  of  3D  motion  once  n  is  known, 
and  (iii)  how  to  compare  the  degree  of  3D  motion  from  one  imaging  interval  to  another. 
These  three  issues  are  discussed  in  the  following  subsections. 


4.1.  Phase  Estimation  Using  Adaptive  Joint  Time-Frequency  Projection 

After  the  translational  motion  compensation,  the  radar  signal  contains  only 
rotational  motion.  To  estimate  the  phase  of  a  prominent  point  scatterer,  we  utilize  the 
adaptive  joint  time-frequency  (AJTF)  projection  technique  discussed  in  [6].  We  begin 
with  the  radar  data  in  the  (range)-(dwell  time)  domain.  Within  a  fixed  range  cell,  the  data 
can  be  written  as 
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(14) 


,=i 

where /c  is  the  center  frequency.  Among  the  Nr  point  scatterers  within  the  range  cell,  we 
express  the  phase  behavior  of  the  strongest  one  as  a  polynomial  function: 

=  +  +  (15a) 

and  consider 

h{t)  =  (to)]  (15b) 

as  a  basis  for  the  radar  signal.  The  phase  parameters  are  then  found  by  searching  for  the 
maximum  projection  from  the  radar  signal  onto  the  basis  function: 

<  /, ,  /2 ,  /j ....  >=  arg  max  \\E(to)h\to)dto\  (16) 

Equation  (16)  means  that  the  phase  function  parameters  are  estimated  to  give  a  maximum 
projection  from  the  radar  data  onto  the  basis  function  for  that  prominent  point  scatterer. 
In  the  search  procedure,  the  first  term  fi  can  be  obtained  by  using  the  fast  Fourier 
transform,  while  all  other  higher  order  terms  fz,  fs,  •••  are  obtained  using  exhaustive 
search.  Figure  5  illustrates  the  process  of  AJTF  phase  estimation.  Figure  5(a)  shows  the 
radar  signal  in  one  range  cell  with  three  point  scatterers  in  the  joint  (dwell  time)-(Doppler 
frequency)  plane.  The  tilted  trajectory  of  the  prominent  point  scatterer  1  implies  there 
exist  higher-order  terms  in  the  phase  function.  Figure  5(b)  shows  the  basis  function  h(tD)- 
During  the  search,  we  change  the  position  (/}),  tilting  (fz)  and  curvature  (fj,  ...)  of  h  until 
the  projection  of  h  onto  the  radar  signal  is  maximized. 
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4.2.  Measure  of  Nonlinearity  between  Two  Phase  Functions 

We  notice  that  in  (8),  the  two  phase  functions  are  formulated  with  a  linear 
relationship  plus  a  nonlinear  residual  part.  After  the  two  phase  functions  are  estimated 
using  the  AJTF  technique,  a  least-squares  fitting  can  be  performed  to  generate  the  best-fit 
linear  part.  The  actual  phases  deviates  from  this  linear  relationship.  The  deviation  n  is 
integrated  over  the  dwell  time  to  represent  the  degree  of  phase  nonlinearity  over  the 
imaging  interval  as  follows: 

=  (17) 

The  process  is  illustrated  in  Figure  6.  The  solid  line  is  the  actual  relationship  between  the 
two  phase  functions  Pi  and  Pa-  The  dotted  line  is  the  linear  approximation  of  the 
relationship.  The  area  of  the  shadowed  region  is  Nn- 

In  a  similar  fashion,  we  define  the  degree  of  3D  motion  as  the  deviation  from  a 
linear  relationship  between  Pand  ^over  the  dwell  interval  as  follows: 

M  =  |l  m(ro )  I  (itp  (18) 

Based  on  (13),  we  see  that  M  and  N12  are  directly  related: 

Nn=PnM  (19a) 

where 


Ay,  -I-  PAz, 


(19b) 


Thus  by  finding  the  observable  A^/a,  we  can  obtain  the  degree  of  3D  motion  M  to  within  a 


proportionality  constant. 
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4.3. 3D  Motion  Comparison  among  Different  Imaging  Intervals 

As  indicated  by  (19),  the  phase  nonlinearity  of  two  point  scatterers  N  is 
proportional  to  the  degree  of  3D  motion  M,  so  we  can  use  the  detected  phase  nonlinearity 
as  a  measure  of  3D  motion.  However,  we  notice  that  the  constant  of  proportionality  is 
dependent  on  the  point  scatterer  positions.  A  problem  arises  when  we  need  to  compare 
the  detection  result  from  one  imaging  interval  to  that  from  another  imaging  interval. 
Since  we  cannot  guarantee  that  we  track  the  same  set  of  points  from  frame  to  frame,  the 
proportionality  constant  can  change  from  frame  to  frame,  and  we  cannot  reliably  observe 
M  from  N  across  frames.  To  overcome  this  difficulty,  we  track  more  than  two  point 
scatterers  within  each  frame  and  compute  Nij  for  each  pairing  of  scatterers  i  andj 
Then  we  generate  an  average  value  <Nij>  from  all  the  possible  phase  relationships. 
From  (19),  we  have 

<  Ny  >=<  fiy>M  (20) 

We  postulate  that,  from  a  statistical  point  of  view,  <fiij>  approaches  a  constant  that  is 
independent  of  frames  if  we  average  over  a  sufficient  number  of  point  scatterers.  If  this 
is  true,  <A/y>  should  become  a  good  indicator  of  M. 

We  test  the  effectiveness  of  this  approach  on  the  detection  result  by  simulation. 
We  input  a  set  of  motion  parameters  and  generate  the  phase  functions  based  on  the  3D 
motion  model.  20  point  scatterers  from  an  airplane  model  is  used.  We  then  randomly 
choose  a  number  of  point  scatterers  and  use  their  phase  functions  to  compute  <Nij>.  We 
examine  how  the  results  vary  as  different  number  of  point  scatterers  is  used.  We  find  that 
the  results  begin  to  converge  after  about  5  scatterers.  Fig.  7  shows  a  plot  of  <Nij>  versus 
the  frame  number  if  we  use  5  point  scatterers  (10  phase  pairs).  If  we  increase  the  number 
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of  point  scatterers  to  10  (45  phase  pairs),  there  is  only  minor  change  in  the  detection 
output.  Therefore,  <Nij>  can  be  used  to  indicate  the  degree  of  3D  motion  given  a 
sufficient  number  of  point  scatterers. 

5.  Results 

To  demonstrate  the  effectiveness  of  the  3D  motion  detection  algorithm,  we  test 
our  algorithm  on  radar  data  from  two  targets.  The  first  target  is  an  aircraft,  which  flew  in 
a  large  clockwise  circle  during  a  9-minute  interval.  We  also  have  access  to  the  target 
motion  data  through  the  GPS  (global  positioning  system)  and  INS  (inertial  navigation 
system)  sensors  carried  on-board  the  aircraft  [14-16].  Figure  8  shows  our  processing  flow 
chart.  The  GPS/INS  data  is  used  to  establish  the  tmth  target  motion.  The  raw  radar  data  is 
used  as  input  to  the  3D  motion  detection  algorithm.  We  can  also  generate  the  ISAR 
images  using  our  AJTF  motion  compensation  algorithm  [6].  We  are  therefore  able  to  both 
compare  the  detection  result  with  the  truth  motion,  and  observe  the  effect  of  the  3D 
motion  on  the  ISAR  image  quality. 

We  first  test  the  3D  motion  detection  algorithm  on  simulated  radar  data.  To 
generate  the  simulation  data,  we  use  the  actual  motion  data  from  the  GPS/INS  sensors  in 
conjunction  with  a  point  scatterer  model.  From  the  aircraft  model,  60  point  scatterers  are 
selected  to  simulate  the  radar  data  based  on  the  actual  motion  data  and  equation  (2).  Five 
range  cells  are  then  chosen  for  phase  analysis  in  the  detection  procedure.  Figure  9(a) 
shows  the  detected  degree  of  3D  motion  for  20  image  frames  from  the  simulated  radar 
data.  For  comparison.  Figure  9(b)  shows  the  degree  of  3D  motion  obtained  based  on  the 
tmth  motion  data.  The  frames  with  significant  3D  motion  are  highlighted  with  circles  and 
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the  frames  with  2D  motion  are  highlighted  with  diamonds.  It  is  seen  that  the  two  results 
agree  fairly  well. 

Next,  we  test  the  detection  algorithm  using  the  actual  radar  measurement  data. 
Figure  10(a)  shows  the  detected  3D  motion  from  the  radar  data  over  20  frames.  The 
corresponding  imaging  interval  for  each  frame  is  2.3  seconds  while  the  total  flight 
duration  is  5  minutes.  The  four  frames  with  the  most  significant  3D  motion  based  on  our 
detection  algorithm  are  labeled  as  circles.  They  are  frames  6,  14,  17  and  18.  Figure  10(b) 
shows  the  degree  of  3D  motion  obtained  based  on  the  truth  motion  data.  We  observe  that 
the  truth  motion  data  indeed  contains  a  high  degree  of  3D  motion  at  those  four  frames 
detected  by  our  algorithm. 

To  further  examine  the  quality  of  the  ISAR  images  when  3D  motion  is  present, 
we  generate  images  using  our  motion  compensation  algorithm  in  Figures  11  to  14.  Figure 
1 1(a)  shows  the  plot  of  6  vs.  ^derived  from  the  truth  motion  data  for  frame  18,  which  is 
a  frame  found  to  contain  substantial  3D  motion.  The  actual  motion  is  shown  in  the  solid 
curve  and  the  dashed  line  is  the  best-fit  2D  motion  approximation.  It  is  clear  that  the  solid 
curve  deviates  significantly  from  the  dashed  line  and  the  actual  motion  cannot  be  well 
approximated  with  2D  motion.  Figure  1 1  (b)  shows  the  resulting  image  obtained  after  the 
motion  compensation,  and  is  blurred  in  the  Doppler  dimension  (vertical  axis).  As 
expected,  the  2D  motion  compensation  algorithm  cannot  focus  all  the  points  due  to  the 
3D  target  motion.  Figures  12  (a)  and  12(b)  show  the  same  conclusion  for  frame  14,  which 
is  another  frame  identified  as  having  significant  3D  motion.  In  Figure  13,  we  show  the 
results  for  frame  2,  which  has  very  little  3D  motion.  As  we  can  see  from  Figure  13(a),  the 
actual  motion  can  be  well  approximated  by  a  line  in  the  ^^plot.  The  image  shown 
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Figure  13(b)  is  well  focused.  In  particular,  the  point  scatterers  on  the  target  show  nearly 
equal  range  and  Doppler  extent,  contrary  to  the  previous  two  images.  The  aircraft 
bodyline  is  clearly  recognizable. 

From  Figure  10,  we  notice  that  there  exists  a  discrepancy  in  frame  11,  where  the 
detection  result  does  not  indicate  any  3D  motion  while  the  truth  motion  data  shows  a 
significant  amount  of  3D  motion.  The  truth  motion  is  shown  in  Figure  14(a),  confirming 
the  presence  of  3D  motion.  One  explanation  is  that  those  prominent  points  used  by  the 
detection  algorithm  lie  nearly  on  a  2D  plane  so  that  they  still  can  be  focused.  As  we  have 
discussed  in  Section  3,  the  2D  model  is  applicable  if  either  the  motion  is  2D  or  the  target 
is  of  2D  in  extent.  It  is  likely  that  the  latter  condition  is  met  for  this  frame.  This  is 
confirmed  by  the  image  shown  in  Figure  14(b).  We  see  that  the  image  quality  is  actually 
not  so  bad.  Therefore,  our  detection  algorithm  objectively  reflects  the  quality  of  the 
images  generated  by  the  2D  motion  compensation. 

A  second  data  set  is  used  to  test  our  3D  motion  detection  algorithm.  This  data  set 
consists  of  the  ISAR  data  collected  from  a  small  ship  on  the  ocean.  Because  of  the 
surface  movement  of  the  sea,  the  target  is  believed  to  have  considerable  3D  motion 
during  the  imaging  intervals.  The  3D  motion  detection  result  is  shown  in  Figure  15(a) 
with  the  peaks  corresponding  to  regions  with  3D  motion.  The  total  data  duration  is  20 
seconds  and  the  imaging  dwell  time  is  0.64  second  per  frame.  For  this  data  set,  reliable 
truth  target  motion  is  not  available.  Instead,  we  generate  the  motion  compensated  images 
shown  in  Figures  15(b)  to  15(d)  to  demonstrate  the  effect  of  3D  motion  on  ISAR  image 
quality.  The  image  frame  with  the  largest  detected  3D  motion,  frames  3,  is  shown  in 
Figure  15(b).  It  is  poorly  focused.  Figure  15(c)  shows  the  image  from  frame  14,  which 
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is  the  frame  with  the  second  highest  detected  3D  motion.  The  frame  with  the  smallest  3D 
motion  based  on  our  algorithm,  frame  20,  is  shown  in  Figure  15(d).  It  shows  a  well- 
focused  ISAR  image.  This  test  confirms  the  effectiveness  of  our  algorithm  in  detecting 
good  imaging  intervals  from  those  imaging  intervals  containing  large  3D  motion. 

6.  Conclusions 

In  this  paper,  we  set  out  to  develop  an  algorithm  to  detect  the  presence  of  3D 
target  motion  from  ISAR  data.  Based  on  the  3D  point  scatterer  model,  we  first  examined 
the  effect  of  3D  motion  on  ISAR  imaging.  It  was  shown  that  the  existing  motion 
compensation  algorithms  could  not  properly  focus  targets  exhibiting  3D  motion  during 
the  imaging  interval.  We  then  derived  an  algorithm  to  blindly  detect  the  degree  of  3D 
target  motion  from  raw  radar  data.  It  is  based  on  measuring  the  linearity  of  phases 
between  two  or  more  point  scatterers  on  the  target.  The  phase  estimation  was 
implemented  using  the  adaptive  joint  time-frequency  technique.  Examples  were  provided 
to  demonstrate  the  effectiveness  of  the  3D  motion  detection  algorithm  with  both 
simulation  and  real  ISAR  data.  The  detection  results  were  corroborated  with  the  truth 
motion  data  from  on-board  motion  sensors  and  correlated  with  the  resulting  ISAR 
images.  With  the  detection  algorithm,  we  have  the  ability  to  distinguish  the  time  intervals 
when  the  target  undergoes  smooth  2D  motion  from  those  containing  more  chaotic  3D 
motion.  As  a  result,  the  good  imaging  intervals  where  focused  images  are  more  easily 
formed  can  be  automatically  selected. 
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(a)  (b) 


Fig.  3  Phase  linearization  achieved  by  rotational  motion  compensation  (a)  Phases  of 

two  point  scatterer  with  2D  rotational  motion,  (b)  Both  phases  are  linearized 
with  rotational  motion  compensation,  (c)  Phases  of  two  point  scatterers  with  3D 
motion,  (d)  Only  one  phase  is  linearized  with  rotational  motion  compensation. 


Fig.  4  Problem  with  a  typical  motion  compensation  algorithm  (a)  Target  undergoes  a  2D 
motion,  (b)  Image  after  translational  motion  compensation,  (c)  Image  after 
rotational  motion  compensation,  (d)  Target  undergoes  a  3D  motion,  (e)  Image  after 
translational  motion  compensation,  (f)  Image  after  rotational  compensation. 
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Fig.5  (a)  (Dwell  time)-(Doppler  frequency  representation  of  radar  signal  in  a 

range  cell  with  three  point  scatterers.  (b)  The  basis  function  that  is  best  matched 
to  the  dominant  point  scatterer  is  found  by  the  AJTF  project  method. 


Fig.  6  Measure  of  the  nonlinearity  of  two  phase  functions. 


Degree  „ 
of  3D 
motion  „ 

0,4 

0.2 

0 


Fig.  7  Effect  of  the  number  of  point  scatterers  used  on  3D 
motion  detection  result. 


Fig,  8  Data  processing  flow  chart. 
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Fig.  10  (a)  Detected  3D  motion  from  aircraft  radar  data, 
(b)  Degree  of  3D  motion  from  truth  motion  data. 
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Fig.  1 1  3D  motion  and  resulting  ISAR  image  (frame  18). 
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Fig.  12  3D  motion  and  resulting  ISAR  image  (frame  14). 
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Fig.  13  2D  motion  and  resulting  ISAR  image  (frame  2) 
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Abstract  In  this  paper,  genetic  algorithms  (GA)  are  proposed  for  inverse  synthetic 
aperture  radar  (ISAR)  imaging  of  an  unknown  target.  The  adaptive  joint  time-frequency 
(AJTF)  analysis  is  used  to  achieve  translation  and  rotational  motion  compensation 
through  motion  parameterization.  GA  is  used  as  an  alternative  to  exhaustive  search  in  the 
parameter  search  process.  While  maintaining  the  same  accuracy,  GA  has  lower 
computational  complexity,  especially  for  targets  with  highly  irregular  motions. 
Simulation  and  measurement  data  results  demonstrate  the  feasibility  of  GA  in  ISAR 
imaging  of  targets  with  high-order  motions. 
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1.  Introduction 

Inverse  synthetic  aperture  radar  (ISAR)  system  usually  collects  radar  data  from  a 
target  moving  on  the  ground,  in  the  air,  or  over  the  ocean.  In  the  ISAR  problem  shown  in 
Figure  1 ,  the  radar  is  stationary  while  the  target  moves  with  both  translation  motion  and 
rotational  motion.  In  the  microwave  frequency  range,  ISAR  has  been  identified  as  an 
effective  tool  for  target  identification  [1]. 
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Target  motion  is  an  essential  part  in  ISAR  imaging.  On  the  one  hand,  target 
motion  is  needed  to  generate  Doppler  (or  cross-range)  resolution  [2].  On  the  other  hand, 
unwanted  motion  causes  image  blurring.  When  the  target  has  uniform  rotational  motion 
only  and  the  radar  data  is  collected  over  a  small  angular  aperture,  a  simple  Fourier 
transform  will  bring  the  raw  radar  data  into  a  two-dimensional  ISAR  image.  However, 
actual  targets  observed  by  operational  radar  rarely  have  such  an  ideal  motion.  Therefore, 
motion  compensation  is  needed  to  generate  focused  ISAR  imagery.  There  exist  many 
different  motion  compensation  algorithms  to  deal  with  target  motion  [3-5].  Most  of  the 
algorithms  start  with  coarse  range  alignment  based  on  the  correlation  of  the  range 
profiles.  Then  the  phase  information  within  one  range  cell  is  utilized  to  achieve  fine 
motion  compensation. 

Phase  estimation  is  critical  in  fine  motion  compensation.  Compared  to  the 
amplitude,  the  phase  of  the  radar  signal  is  much  more  sensitive  to  the  change  in  range. 
Based  on  the  concept  of  signal  parameterization  reported  in  [6,  7],  an  adaptive  joint  time- 
frequency  (AJTF)  algorithm  was  proposed  in  [8]  for  phase  estimation  of  a  prominent 
point  scatterer.  In  this  method,  the  target  motion  is  modeled  as  a  polynomial  function  and 
an  exhaustive  search  procedure  is  used  to  find  the  motion  parameters  that  are  embedded 
in  the  phase  of  the  prominent  point  scatterer.  While  the  performance  of  this  algorithm  is 
very  good  [9],  the  main  bottleneck  in  this  procedure  is  the  computational  complexity 
associated  with  the  parameter  search.  When  the  target  motion  is  highly  irregular,  the 
number  of  parameters  needed  to  model  the  motion  becomes  large  and  the  use  of  the 
exhaustive  search  becomes  prohibitively  expensive. 
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In  this  paper,  our  objective  is  to  reduce  the  computation  time  associated  with  the 
motion  parameter  search  in  the  AJTF  procedure.  Our  proposed  approach  is  to  incorporate 
genetic  algorithms  (GA)  [10]  into  the  AJTF  search  process.  (Some  preliminary  work  on 
this  topic  was  reported  in  [11].)  In  contrast  with  conventional  optimization  methods,  GA 
is  a  population-based,  statistical  search  technique.  It  borrows  such  concepts  as  inheritance 
and  mutation  from  the  biological  evolution  process  [12].  As  a  global  optimization 
technique,  GA  is  known  to  be  very  easy  to  implement  and  applicable  to  many  design  and 
inverse  problems  [13]. 

This  paper  is  organized  as  follows.  In  Sections  2  and  3,  we  outline  the 
methodology.  The  AJTF  analysis  for  ISAR  motion  compensation  is  described  in  Section 
2.  Genetic  algorithms  are  introduced  in  Section  3.  The  next  two  sections  include  results 
and  analysis.  In  Section  4,  simulations  with  point  scatterers  are  provided  to  validate  the 
use  of  GA  for  phase  estimation.  Measurement  data  processing  results  are  shown  in 
Section  5.  Finally,  conclusions  are  given  in  Section  6. 

2.  ISAR  Motion  Compensation  Using  the  Joint  Time-Frequency  Projection 

A  two-dimensional  point  scatterer  model  relates  the  radar  data  to  a  moving  target 
through  the  following  equation 

=  )] }  (1) 

where  /  is  the  frequency  and  to  is  the  dwell  time.  In  this  model,  the  radar  data  is 
comprised  of  the  sum  of  responses  from  a  collection  of  point  scatterers.  (xi,yi)  denotes  the 
point  scatterer  position  while  Oj  denotes  the  scatterer  strength.  The  target  motion  includes 
both  translation  motion  rfto)  and  rotational  motion  ^to)- 
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After  range  compression  is  done  by  Fourier  transforming  the  data  with  respect  to 
frequency,  the  original  radar  data  is  converted  into  a  set  of  range  profiles.  Range 
alignment  is  then  carried  out  via  amplitude  correlation  of  the  range  profiles  to  place  all 
the  point  scatterers  in  their  correct  range  bins.  The  radar  signal  through  one  range  cell  r 
can  thus  be  expressed  in  the  form  of 

^r(^D)  =  exp{-7-^^[Ar(fD)  +  x,.  +  (2) 

/  c 

where  fo  is  the  center  frequency  of  the  radar  and  the  index  includes  only  those  point 
scatterers  in  the  particular  range  cell.  The  residual  translation  motion  is  depicted  as 
Ar(tD)-  After  such  coarse  alignment  procedure,  the  residual  translation  motion  is  smaller 
than  the  range  resolution.  However,  it  can  still  be  larger  than  the  radar  wavelength. 
Therefore,  it  is  important  to  include  this  term  in  (2). 

Both  of  the  residual  translation  motion  Ar(to)  and  the  rotational  motion  (pita)  can 
be  expanded  into  polynomial  functions  of  the  dwell  time  as 

where  any  coefficients  beyond  the  first  order  are  detrimental  to  ISAR  image  formation. 
To  solve  the  ISAR  motion  compensation  problem,  we  need  to  determine  these  motion 
parameters  and  to  remove  the  unwanted  nonlinear  phase  terms  from  the  radar  data. 

This  task  can  be  accomplished  using  the  AJTF  procedure  [8].  The  essential  idea 
of  this  procedure  is  to  find  the  basis  function  that  most  resemble  the  strongest  signal 
component  in  equation  (2).  For  our  problem,  a  basis  function  in  the  form  of 

h(t^ )  =  exp[-7  (/,to  +  f^tl  + 
c 
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is  used.  The  best  basis  is  found  by  searching  for  parameters //,/2,/i,  ...  that  maximize  the 
projection  from  the  radar  signal  onto  the  basis,  i.e., 

<  fv  fi^  /i.-"  >=  argmax  k  E/toX  >' 

where  the  projection  is  formulated  as  the  inner  product  of  the  two  functions  as 

<  EMMto)  >=  lEMh\tD)dto  (6) 

In  (5),  the  linear  coefficient//  can  be  found  efficiently  with  the  fast  Fourier  transform 
(FFT).  Coefficients  for  the  nonlinear  phase  terms, /2,/i,  ...,  must  be  determined  through  a 
more  time-consuming  search. 

Figure  2  depicts  the  concept  of  AJTF  processing  for  ISAR  translation  motion 
compensation.  After  range  alignment,  the  range  profiles  are  shown  in  Figure  2a.  Range 
cell  i  contains  a  prominent  point  scatterer  and  Figure  2b  shows  the  joint  (dwell  time)- 
(Doppler  frequency)  representation  of  the  radar  signal  through  this  range  cell.  A  signal 
with  linear  phase  behavior  shows  up  as  a  straight  horizontal  line  in  the  joint  time- 
frequency  plane.  Both  scatterers  shown  in  Figure  2b  are  tilted  and  curved  due  to  the 
higher-order  residual  motion.  Also  shown  in  the  figure  is  the  basis  function  h.  To  get  the 
maximum  projection,  different  values  of  fhfi,  fs,  •••  are  tried  until  the  basis  function 
approaches  the  strongest  scatterer.  After  the  motion  parameters  of  the  prominent  point 
scatterer  are  found  in  this  manner,  we  can  carry  out  the  translation  motion  compensation 
by  multiplying  the  radar  data  with  the  conjugate  of  this  basis.  Since  all  the  point 
scatterers  share  the  same  translation  motion  in  equation  (2),  this  operation  will  remove 
the  translation  motion  of  the  whole  target. 

After  the  translation  motion  compensation,  we  can  rewrite  equation  (2)  as 
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(7) 


^r(^D)  =  exp{-7^^[Ax,.  +  +b/o...)]} 

i  c 

where  (Axi,Ayi)  is  the  differential  position  of  the  scatterer  relative  to  the  prominent 
point  scatterer  chosen  for  translation  motion  compensation  in  the  previous  step.  To 
accomplish  the  rotational  motion  compensation,  a  second  prominent  point  scatterer, 
usually  in  a  range  cell  different  from  the  first  one,  is  chosen  for  phase  analysis.  The  same 
search  procedure  as  described  before  is  used  to  find  the  rotational  motion  parameters  in 
(7).  Then  to  achieve  rotational  motion  compensation,  we  can  simply  define  a  new  dwell 
time  variable  t'o  '■ 

4  =  h  +  +  f’ifo  +  •••) 

and  resample  the  radar  data  uniformly  in  terms  of  this  new  time  variable: 

^r(4 )  =  Z  exp[-y^^(Ax,.  +  )]  (9) 

i  c 

This  resampling  procedure  is  applied  to  the  entire  data  set  to  linearize  the  phase  of  the  all 
the  point  scatterers  on  the  target.  From  equation  (9),  it  is  clear  that  Fourier  transforming 
the  data  with  respect  to  t'p  will  resolve  the  point  scatterers  in  the  cross  range  dimension. 

3.  Use  of  Genetic  Algorithms  for  Phase  Parameter  Search 

As  we  have  discussed  in  the  last  section,  ISAR  motion  compensation  can  be 
formulated  as  a  parameterization  process  for  both  translation  motion  and  rotational 
motion.  In  [8],  a  brute-force  search  procedure  is  employed  to  carry  out  the 
parameterization.  This  means  that  we  exhaustively  search  the  solution  space  for  the 
maximum  projection.  Genetic  algorithms  (GA)  are  investigated  here  to  search  for  the 
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motion  parameters  to  reduce  the  computation  time.  We  should  point  out  that  although  a 
structured  tree-search  is  an  easy  and  straightforward  way  to  decrease  the  computational 
complexity,  it  does  not  always  lead  to  a  global  optimum. 

GA  is  a  global  optimization  method  based  on  concepts  from  ecological  systems 
[10,  12,  14].  The  flowchart  of  a  typical  GA  process  is  illustrated  in  Figure  3.  It  starts  by 
setting  up  the  parameters  for  both  the  physical  problem  and  the  GA  implementation.  GA 
operates  on  a  population  of  many  individuals.  The  initial  population  is  randomly 
generated  within  the  searching  space.  The  goodness  of  the  solution  is  then  evaluated 
based  on  an  objective  function.  For  our  problem,  the  objective  function  is  defined  as  the 
projection  magnitude  from  the  radar  signal  onto  the  basis  function  (see  equation  (5)).  If 
we  are  satisfied  with  the  solution,  the  process  is  done.  Otherwise,  a  new  generation  is 
produced  for  the  next  evaluation-reproduction  iteration.  To  make  up  individuals  in  the 
new  generation,  some  good  parents  are  selected  from  the  previous  generation  and  two 
operations  called  crossover  and  mutation  are  applied  to  produce  children  in  the  next 
generation.  Whether  or  not  crossover  and  mutation  occur  is  determined  randomly.  The 
crossover  and  mutation  probabilities  are  chosen  based  on  the  tradeoff  between  two 
conflicting  requirements.  Increasing  the  variation  in  the  new  generation  brings  a  chance 
for  better  solutions,  but  it  tends  to  lose  the  features  of  the  good  solutions  from  the 
previous  generation. 

Roughly  speaking,  there  are  two  kinds  of  GA.  One  is  binary-coded  GA  [10].  The 
other  is  real-coded  GA  [15].  In  the  former,  the  physical  parameters  to  be  searched  are 
first  discretized  into  binary  bits.  There  is  a  one-to-one  mapping  between  a  physical 
parameter  C  and  its  iV-bit  binary  representation  c/,  cj, ...,  cn  as  follows: 
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where  the  [Cmin,  Cmax]  is  the  valid  search  space  for  C.  A  candidate  solution  of  the  problem 
is  expressed  in  the  form  of  a  chromosome,  which  is  the  collection  of  bits  representing  all 
the  parameters.  For  crossover,  a  random  break  point  in  the  chromosome  is  picked.  The 
bits  before  the  point  from  one  parent  are  combined  with  the  bits  after  the  point  from 
another  parent  to  form  one  child.  Another  child  is  generated  in  the  reverse  fashion.  For 
mutation,  a  single  bit  is  randomly  picked  and  its  valued  is  inverted.  The  crossover  and 
mutation  operations  for  binary-coded  GA  are  depicted  in  Figure  4a. 

In  real-coded  GA,  there  is  no  coding-decoding  process  and  the  algorithm  directly 
operates  on  the  physical  parameters.  For  crossover,  a  linear  combination  is  usually  used 


as  follows 


J  C,  =  «P,  +  (1  -  a)P2 

[Cj  =  (1  -  a)Pi  +  «p2 


where  two  children  (C/,  C2)  are  reproduced  from  two  parents  {PhPi)-  a  is  a  random 
number  between  0  and  1  to  ensure  that  the  new  parameters  will  not  fall  out  of  range.  For 
mutation,  a  child  C  that  is  different  from  the  parent  P  is  needed.  For  this  purpose,  a 
solution  Pe  is  picked  up  randomly  from  the  searching  space  and  linearly  combined  with  P 
to  generate  C  as  in  equation  (11).  The  crossover  and  mutation  operations  for  real -coded 
GA  are  depicted  in  Figure  4b. 

The  basic  theory  in  GA,  the  schemata  theory,  seems  to  favor  the  use  of  binary- 
coded  GA  [14].  Most  work  on  GA  has  followed  this  path.  Recently,  researchers  have  also 
experimented  with  real-coded  GA  and  have  observed  some  advantages  in  convenience 


8 


and  accuracy  [16].  In  the  next  section,  we  will  test  both  binary-coded  and  real-coded  GA 
in  our  phase  estimation  problem. 

The  GA  process  is  usually  stopped  using  criteria  based  on  the  performance  of  the 
available  solutions  in  the  present  generation.  In  our  case,  we  do  not  know  the  true 
maximum  projection  value.  Therefore,  we  choose  to  stop  the  GA  process  when  the 
projection  value  does  not  increase  after  a  certain  number  of  generations. 


4.  Point  Scatterer  Simulation 

Point  scatterer  simulations  are  first  used  to  test  the  use  of  GA  for  ISAR  motion 
compensation.  In  the  first  example,  we  assume  two  point  scatterers  with  amplitudes  1  and 
0.2.  They  are  located  within  one  range  cell  and  contain  third-order  translation  motion 
(i.e.,  the  coefficients  a/,  <22  and  as  in  equation  (3)  are  significant  while  all  higher-order 
coefficients  are  zero).  We  run  both  binary-coded  GA  and  real-coded  GA  to  search  for  <22 
and  as  for  this  simple  phase  estimation  problem.  The  population  size  is  50.  In  both  cases, 
the  crossover  probability  is  0.8  and  the  mutation  probability  is  0.3.  For  exhaustive  search, 
the  search  for  02  and  as  is  carried  out  on  a  128  by  128  grid.  We  choose  a  7-bit 
representation  in  binary-coded  GA  and  search  in  the  same  discrete  space  as  in  the 
exhaustive  search.  The  actual  objective  function  surface  is  plotted  in  Figure  5a.  We 
observe  many  local  maxima,  indicating  this  would  be  a  challenging  problem  for  a  local 
optimization  method.  Figure  5b  shows  the  GA  convergence  curve,  with  the  real-coded 
GA  result  in  solid  line  and  the  binary-coded  GA  result  in  dashed  line.  Since  GA  is  a 
statistical  approach,  we  do  not  get  exactly  the  same  result  from  each  GA  run.  To  decrease 
the  statistical  variation,  the  results  in  Figure  5b  are  obtained  by  averaging  over  20  GA 
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runs.  We  can  see  that  after  about  150  generations  the  two  projection  curves  nearly 
converge  to  the  truth  value  of  1 .  We  also  observe  that  real-coded  GA  produces  a  slightly 
higher  projection  value.  Figure  5c  shows  the  resulting  phase  from  a  single  GA  run  after 
200  generations.  The  estimated  phase  from  binary-coded  GA  is  plotted  in  circles,  the 
phase  from  real-coded  GA  in  crosses,  and  the  original  phase  function  in  solid  line.  We 
see  very  good  agreement  among  the  three  results,  meaning  good  accuracy  from  the  two 
GA  results. 

In  the  second  example  we  compare  the  computational  complexity  of  GA  to 
exhaustive  search  for  different  orders  of  motion.  As  we  have  pointed  out  earlier,  the  main 
problem  with  exhaustive  search  for  motion  parameter  extraction  is  the  computational 
load.  This  problem  becomes  more  acute  when  the  order  of  the  motion  is  high.  Again,  we 
use  two  point  scatterers  with  amplitudes  1  and  0.2.  We  generate  the  radar  data  from  these 
two  point  scatterers  with  some  preset  motion.  We  then  apply  exhaustive  search,  binary- 
coded  GA,  and  real-coded  GA  for  the  phase  estimation  problem  with  different  orders  of 
motion.  The  same  GA  parameters  as  in  the  previous  example  are  used  and  the  results  are 
averaged  over  20  runs.  The  exhaustive  search  is  known  to  have  an  exponential 
complexity  of  0(exp(n)).  As  expected,  the  resulting  computation  time  in  logarithm  scale 
shows  up  as  a  straight  line  in  Figure  6a.  For  GA,  no  theoretical  complexity  formulation  is 
available  in  general.  (The  complexity  of  0(n  log  n)  is  claimed  for  selected  test  functions 
in  [17].)  We  run  both  binary-coded  GA  and  real-coded  GA  up  to  sixth-order  motion  (i.e., 
search  for  5  parameters)  and  plot  the  results  in  Figure  6a.  It  is  observed  that  both  binary- 
coded  GA  and  real-coded  GA  have  much  lower  complexity  than  exhaustive  search.  The 
difference  in  complexity  between  the  two  GA  is  only  minor.  The  projection  values  from 
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binary-coded  GA  and  real-coded  GA  are  next  plotted  in  Figure  6b  as  circles  and  crosses, 
respectively.  We  see  that  the  real-coded  GA  results  are  closer  to  the  truth  value  of  1.0 
than  the  binary-coded  GA  results,  especially  for  higher-order  motions.  Since  binary- 
coded  GA  searches  on  a  finite  grid  as  in  exhaustive  search,  it  can  never  get  solutions  that 
surpass  the  exhaustive  search  result.  On  the  other  hand,  real-coded  GA  has  the  ability  to 
search  for  any  real  values  within  the  search  range.  Consequently,  real-coded  GA  has 
more  chance  of  finding  a  better  solution.  The  same  trend  is  also  observed  with 
measurement  data  and  will  be  discussed  further  in  the  next  section. 

In  the  last  simulation  example,  we  show  that  the  correct  phase  estimation  from 
GA  can  be  utilized  to  achieve  ISAR  motion  compensation.  For  this  purpose,  four  point 
scatterers,  with  two  in  one  range  cell  and  two  in  another  range  cell,  are  assumed  to  carry 
both  translation  and  rotational  motion  up  to  third  order.  The  aligned  range  profiles  are 
shown  in  Figure  7a.  The  dwell  time  vs.  Doppler  frequency  spectrograms  of  the  original 
radar  signal  in  range  cells  1  and  2  are  shown  in  Figures  7b  and  7c  respectively.  Figure  7d 
shows  the  resulting  cross  range  image  using  the  Fourier  transform  directly.  Without 
motion  compensation,  the  four  point  scatterers  are  not  well  resolved.  After  we  apply  real- 
coded  GA  to  the  data  in  range  cell  1,  the  estimated  phase  function  is  shown  in  Figure  8a. 
This  corresponds  to  the  translation  motion  of  the  target.  After  removing  this  motion,  the 
resulting  signal  has  spectrograms  for  range  cells  1  and  2  shown  in  Figures  8b  and  8c 
respectively.  As  expected,  only  one  point  is  focused  in  Figure  8d.  We  next  estimate  the 
phase  of  a  second  strong  point  scatterer  by  mnning  GA  for  range  cell  2.  The  resulting 
phase  function  is  shown  in  Figure  9a.  Time  resampling  is  done  to  linearize  this  curve. 
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Now  all  the  point  scatterers  appear  as  horizontal  lines  in  Figures  9b  and  9c.  Figure  9d 
shows  the  resulting  image.  The  four  point  scatterers  are  all  well  resolved  in  this  figure. 

5.  Measurement  Data  Processing 

We  next  apply  GA  on  some  measurement  data.  The  data  were  collected  from  an 
in-flight  aircraft  using  a  ground  radar.  128  pulses  are  processed  to  form  an  ISAR  image. 
This  corresponds  to  an  imaging  interval  of  about  2.5  seconds.  GA  is  evaluated  for  fine 
motion  compensation.  For  the  first  data  set,  the  image  without  any  fine  motion 
compensation  is  shown  in  Figure  10a.  This  image  is  unfocused  due  to  the  residual 
motion.  We  select  a  range  cell  with  a  dominant  scatterer  (range  cell  79)  and  apply  GA  to 
estimate  the  phase  based  on  a  third-order  translation  motion  model.  The  resulting  images 
from  binary-coded  and  real -coded  GA  are  shown  in  Figures  10b  and  10c,  respectively. 
We  observe  that  the  two  GA  images  are  as  focused  as  the  reference  image  in  Figure  lOd 
obtained  using  exhaustive  search.  The  corresponding  projection  values  are  2401  and 
2595,  as  compared  to  the  exhaustive  search  result  of  2401.  We  continue  this  comparison 
using  19  other  imaging  intervals.  Figure  11a  shows  the  projection  values  (normalized 
with  respect  to  the  exhaustive  search  result)  for  the  20  frames.  For  19  out  20  frames,  real- 
coded  GA  gives  larger  projection  values  than  exhaustive  search.  The  resulting  images  are 
either  on  par  or  slightly  better  focused  than  those  obtained  from  exhaustive  search.  For 
binary-coded  GA,  7  frames  have  lower  projection  values  and  are  of  inferior  image  quality 
to  those  from  exhaustive  search.  The  other  13  frames  have  exactly  the  same  projection 
values  as  the  exhaustive  search  result.  Similar  to  our  conclusion  earlier  based  on  the 
simulation  data,  our  experience  with  the  measurement  data  indicates  that  real-coded  GA 
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consistently  outperforms  binary-coded  GA  in  terms  of  accuracy.  The  computation  time 
using  Matlab  codes  on  a  Pentium  III  750MHz  PC  is  summarized  in  Figure  1  lb.  While 
there  is  little  change  in  the  computational  time  for  exhaustive  search  from  one  frame  to 
another,  the  times  for  binary-coded  and  real-coded  GA  exhibit  large  variations  in  these 
single  run  results  due  to  the  statistical  nature  of  GA.  The  average  times  for  the  binary- 
coded  and  real-coded  GA  are  19.5  seconds  and  11.5  seconds,  respectively.  This  is 
compared  to  the  45.5  seconds  from  exhaustive  search.  Finally,  we  note  that  for  these  20 
frames  the  third-order  model  is  adequate  to  model  the  translation  motion.  Inclusion  of 
higher-order  translation  motion  or  rotational  motion  does  not  improve  the  image  quality 
for  this  data  set. 

For  a  second  data  set,  we  first  apply  third-order  translation  motion  compensation 
using  real-coded  GA.  The  resulting  image  is  shown  in  Figure  12a.  It  is  seen  that  the 
selected  dominant  point  scatterer  at  range  cell  64  is  not  well  focused.  This  means  that  the 
target  contains  higher  motion  that  cannot  be  fully  compensated  by  the  third  order  motion 
model.  This  fact  is  also  revealed  by  the  spectrogram  of  the  compensated  signal  in  Figure 
12b,  as  we  observe  a  curved  JTF  trajectory  for  the  scatterer.  Next,  a  fourth-order  motion 
model  is  tried  and  the  real-coded  GA  result  is  shown  in  Figure  13a.  From  this  figure,  the 
reference  point  scatterer  is  better  focused  and  the  spectrogram  of  the  signal  is 
straightened  in  Figure  13b.  Fifth-order  translation  motion  is  also  attempted,  and  the  result 
does  not  show  much  improvement.  However,  we  observe  that  other  point  scatterers  in 
Figure  13a  are  still  smeared.  This  is  confirmed  by  the  spectrogram  of  the  signal  at  a 
different  range  cell  (number  71)  shown  in  Figure  13c,  which  shows  that  the  signal 
trajectory  is  curved  in  the  spectrogram.  Thus  rotational  motion  must  exist  in  this  data.  We 
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next  apply  fourth-order  rotational  motion  compensation  using  real-coded  GA.  As  shown 
in  Figure  14a,  the  whole  target  is  much  better  focused  after  the  compensation.  The 
spectrograms  of  the  signal  at  both  range  cells  become  straightened  in  Figures  14b  and 
14c.  While  real-coded  GA  takes  45  seconds  for  the  phase  estimation  problem,  the 
computation  time  for  fourth-order  exhaustive  search  is  estimated  to  be  over  50  minutes 
based  on  the  complexity  curve  in  Figure  6a.  Therefore,  the  time  savings  of  GA  over 
exhaustive  search  is  quite  significant  in  this  real-world  example. 

6.  Conclusions 

In  this  paper,  genetic  algorithms  have  been  applied  to  ISAR  motion 
compensation.  Based  on  the  adaptive  joint  time-frequency  analysis,  GA  is  used  in  the 
phase  estimation  of  prominent  point  scatterers  on  the  target.  The  resulting  parameterized 
phases  are  then  used  for  translation  and  rotational  motion  compensation.  Both  binary- 
coded  GA  and  real-coded  GA  have  been  implemented  and  tested  using  simulation  and 
measurement  data.  It  is  found  that  real-coded  GA  outperforms  binary-coded  GA  in  terms 
of  accuracy  in  the  phase  estimation  problem.  It  is  also  shown  that  the  computational 
complexity  of  the  GA  search  is  much  less  than  that  of  exhaustive  search.  The  time 
savings  can  become  especially  significant  when  the  target  exhibits  highly  irregular 
motion. 
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Fig.  1 .  Geometry  of  the  IS  AR  imaging  problem. 
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(b) 


Fig.  2.  The  concept  of  AJTF  analysis  for  ISAR  motion  compensation,  (a)  Correlator- 
aligned  range  profiles  with  a  strong  point  scatterer  in  range  cell  i.  (b)  Joint 
(dwell  time)-(Doppler  frequency)  representation  of  the  signal  in  range  cell  i. 
The  basis  function  h  that  is  best  matched  to  the  strongest  scatterer  is 
determined  via  parameter  search. 
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Fig.  3.  Flowchart  of  GA. 
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Fig.  4.  Examples  of  crossover  and  mutation  operations  in  binary  and  real-coded 
GA.  (a)  Binary-coded  GA.  (b)  Real-coded  GA. 
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Fig.  5.  GA  for  third-order  phase  estimation,  (a)  Multi-modal  objective  function,  (b) 
GA  convergence  curves,  (c)  GA-estimated  phase  compared  to  the  original 
truth  phase. 
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Fig.  6.  Performance  of  GA  compared  to  exhaustive  search,  (a) 
Computational  complexity  as  a  function  of  the  number  of 
parameters,  (b)  Accuracy. 


6 


pulse 

No. 


Doppler 

frequency 


dwell  time 


Fig.  7.  Point  scatterer  testing  using  four  scatterers  in  two  range  cells,  (a)  Range 
profiles  prior  to  the  fine  motion  compensation,  (b)  Spectrogram  of  the  signal 
in  range  cell  1 .  (c)  Spectrogram  of  the  signal  in  range  cell  2.  (d)  Resulting 
image  via  FFT  processing. 


Fig.  8.  Translation  motion  compensation  using  real-coded  GA.  (a)  GA-estimated  phase 
due  to  the  translation  motion,  (b)  Spectrogram  of  the  signal  in  range  cell  1.  (c) 
Spectrogram  of  the  signal  in  range  cell  2.  (d)  Resulting  image. 
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Fig.  9.  Rotational  motion  compensation  using  real-coded  GA.  (a)  GA-estimated  phase 
due  to  the  rotational  motion,  (b)  Spectrogram  of  the  signal  in  range  cell  1,  (c) 
Spectrogram  of  the  signal  in  range  cell  2.  (d)  Resulting  image. 
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Fig.  10  Translation  motion  compensation  applied  to  measurement  data  using  a 
third-order  motion  model,  (a)  Image  before  fine  motion  compensation, 
(b)  Binary-coded  GA  result,  (c)  Real-coded  GA  result,  (d)  Exhaustive 
search  result. 
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Fig.  11.  Comparison  of  the  performance  of  GA  and  exhaustive  search  using 
measurement  data,  (a)  Projection  value,  (b)  Computation  time. 
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Fig.  12.  Motion  compensation  example  for  another  measurement  data  set.  (a) 
Image  after  third-order  translation  motion  compensation  using  real- 
coded  GA.  (b)  Spectrogram  of  the  signal  in  range  cell  64. 
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Fig.  13.  Higher-order  translation  motion  compensation,  (a)  Image  after  fourth- 
order  translation  motion  compensation  using  real-coded  GA.  (b) 
Spectrogram  of  the  signal  in  range  cell  64.  (c)  Spectrogram  of  the 
signal  in  range  cell  71. 
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Fig.  14.  Higher-order  rotational  motion  compensation,  (a)  Image  after  fourth-order 
rotational  motion  compensation  using  real-coded  GA.  (b)  Spectrogram  of 
the  signal  in  range  cell  64.  (c)  Spectrogram  of  the  signal  in  range  cell  71. 
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Abstract:  The  genetic  algorithm  is  combined  with  a  moment  method 
solver  to  carry  out  the  shape  inversion  of  two-dimensional  metallic 
objects.  The  binary  bitmap  discretization  is  used  in  conjunction  with 
geometrical  median  filtering  to  describe  the  shape  of  the  object. 
Results  of  the  shape  inversion  using  this  algorithm  are  presented  for 
several  Ipswich  objects  based  on  measurement  data. 

Introduction:  Inverse  synthetic  aperture  radar  (ISAR)  imaging  [1]  is  a 
linearized  form  of  electromagnetic  inverse  scattering.  While  the  algorithm  is 
fast  and  robust  in  obtaining  the  approximate  shape  of  an  object,  it  suffers  from 
resolution  limitation  and  image  artifacts  due  to  multiple  scattering  phenomena. 
Rigorously  solving  the  electromagnetic  inverse  scattering  problem,  on  the  other 
hand,  is  much  more  challenging.  Recently,  some  researchers  have  explored  the 
use  of  genetic  algorithms  (GA)  together  with  computational  electromagnetic 
solvers  to  attack  the  inverse  scattering  problem  [2-6].  In  this  letter,  we  present 
results  of  using  GA  in  conjunction  with  a  method  of  moment  (MoM)  solver  to 
invert  metallic  objects  from  the  Ipswich  measurement  data  set  [7].  In  particular. 


we  focus  our  attention  on  concave  metallic  objects  with  strong  multiple 
scattering  effects. 

In  the  inversion  of  two-dimensional  (2D)  objects,  two  tj^es  of  the 
geometrical  descriptions  have  been  used:  the  Fourier  series  scheme  [2,4]  and 
the  binary  bitmap  discretization  [3,5].  The  Fourier  series  scheme  is  very 
efficient  in  representing  smooth  convex  shapes.  However,  it  does  not  work  well 
for  objects  with  highly  concave  shapes  or  disconnect  parts.  The  binary  bitmap 
discretization  is  a  more  general  way  to  represent  arbitrary  2D  shapes.  Its  main 
drawback  is  the  larger  degrees  of  freedom  required  to  accurately  model  simple 
shapes.  More  recently,  cubic  B-splines  were  also  investigated  as  a  way  to 
accurately  represent  complex  shapes  [6].  In  this  work,  we  use  the  binary 
bitmap  approach  to  discretize  the  search  space.  To  constrain  the  problem,  a 
geometrical  median  filter  is  applied  to  create  more  realizable  shapes.  A  cost 
function  is  defined  as  the  difference  between  the  measured  and  the  computed 
scattered  fields  from  each  assumed  shape.  The  inverse  problem  is  then  cast  into 
a  minimization  problem  and  genetic  algorithm  is  applied  to  minimize  this  cost 
function.  Three  Ipswich  objects,  the  triangular  cylinder,  the  dihedral  and  the 
circular  cavity  Eire  inverted  by  using  this  scheme.  Results  based  on  both  the 
MoM-simulated  field  and  the  measured  data  are  presented. 

GA  Scheme:  In  an  inverse  problem  involving  metallic  objects,  the  measured 
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scattering  field  E  from  the  object  is  known  while  the  shape  and  size  of  the 
object  are  unknown.  We  apply  a  method  of  moment  solution  based  on  the 
electric  field  integral  equation  to  obtain  the  rigorously  computed  scattered  field 


Cdl 

E  from  each  assumed  shape.  To  evaluate  the  performance  of  each  shape,  we 

mea 

define  the  cost  function  as  the  root-mean-square  (RMS)  error  between  E  and 

The  genetic  algorithm  is  applied  as  the  searching  tool  to  minimize  the 
cost  function.  In  our  GA  implementation,  the  initial  generation  is  produced 
randomly  and  each  object  shape  is  encoded  into  an  NxN  binary  array  with  ones 
representing  metal  and  zeros  representing  free  space.  A  2D  median  filter  is  used 
as  a  low-pass  filter  to  eliminate  unrealistic  shapes  consisting  of  isolated  cells. 
With  a  fixed  window  size  of  MxM,  the  median  filter  slides  through  every  cell 
of  the  binary  array  and  sets  the  cell  to  one  if  the  sum  of  the  cell  values  within 
the  window  is  greater  than  or  equal  to  M^/2  and  zero  otherwise.  The  median 
filter  is  applied  once  every  several  generations  so  that  the  isolated  cells  are 
cleared  from  time  to  time,  while  the  new  features  created  by  the  mutation  and 
crossover  operations  have  a  chance  to  survive  in  the  population. 

The  cost  function  for  each  member  is  then  calculated  and  the  shapes 
with  low  cost  values  are  selected  as  parents  to  produce  the  next  generation.  A  2- 
point  crossover  scheme  involving  three  selected  parents  is  used.  The  process 
selects  three  parents  and  divides  each  parent  into  three  parts.  The  three  parent 
shapes  are  then  intermingled  to  produce  three  children  shapes.  This  crossover 
scheme  exhibits  a  more  disruptive  characteristic  for  regeneration  than  the 
conventional  one-point  or  two-point  crossover.  It  serves  to  counteract  against 
the  median  filtering  effect.  The  mutation  operation,  which  is  applied  to  the 
individual  array  cells,  inverts  the  cell  according  to  a  preset  mutation  rate.  The 


selection,  crossover  and  mutation  process  is  iterated  until  the  lowest  eost 
function  in  the  population  reaches  a  sufficiently  small  threshold  or  when  the 
cost  function  does  not  decrease  any  more. 

Results:  The  Ipswich  measurement  data  was  taken  at  a  single  frequency  of 
lOGHz  in  the  bistatic  configuration.  There  were  a  total  of  36  transmitter 
positions  around  the  object  and  18  receiver  locations  for  each  transmitter 
position  (Fig.  1).  Fig.  2a  shows  the  shape  and  size  of  three  metallic  Ipswich 
objects  selected  for  inversion,  namely,  the  triangular  cylinder,  the  dihedral  and 
the  circular  cavity.  They  are  labeled  as  Ips009,  Ips004  and  IpsOl  1,  respectively. 
For  Ips009  and  IpsOl  1,  the  electric  field  is  parallel  to  the  axis  of  the  target.  For 
Ips004,  the  electric  field  is  perpendicular  to  the  axis  of  the  target. 

We  first  tested  the  inversion  algorithm  using  MoM-simulated  field  data 
as  the  input.  In  all  of  the  reconstructions,  the  number  of  the  population  was  set 
to  100  and  the  crossover  and  mutation  rates  were  set  to  0.8  and  0.2, 
respectively.  The  search  area  was  chosen  to  be  15cmxl5cm  for  Ips004  and 
12cmxl2cm  for  Ips009  and  IpsOl  1.  The  number  of  cells  within  this  area  was 
set  to  20x20.  The  reconstructed  results  in  Fig.  2b  show  the  final  inverted 
shapes  of  the  three  objects.  We  observe  that  the  final  shapes  are  in  fairly  good 
agreement  with  the  real  shapes.  The  final  RMS  costs  for  the  three  objects  were 
found  to  be  0.03,  0.38  and  0.18,  respectively.  The  dihedral  and  the  circular 
cavity  contain  strong  multiple  scattering  and  yet  their  inverted  shapes  closely 
resemble  the  correct  objects.  Results  for  these  targets  were  also  generated  using 


the  traditional  imaging  method  and  they  showed  strong  image  artifacts  due  to 
multiple  scattering. 

Next,  we  applied  the  inversion  algorithm  to  the  actual  measured  data. 
Fig.  2c  shows  the  final  reconstruction  shapes.  As  we  can  see,  the  inverted  shape 
is  good  for  the  triangular  cylinder,  which  has  no  multiple  scattering  effects.  For 
the  dihedral,  the  reconstructed  shape  is  not  continuous,  but  is  quite  similar  to 
the  real  object.  The  circular  cavity  shows  the  most  discrepancy  with  the  real 
shape.  The  exterior  and  the  opening  of  the  cavity  are  correctly  inverted,  while 
the  interior  part  of  the  cavity  shape  is  not  as  satisfactory. 

Interestingly,  in  all  three  inversions,  we  found  that  the  final  shape 
obtained  by  GA  has  a  lower  cost  value  than  that  of  the  exact  shape  (0.40  vs. 
0.58  RMS  error  for  the  triangular  cylinder;  0.82  vs.  0.92  for  the  dihedral  and 
0.55  vs.  0.73  for  the  circular  cavity).  This  is  caused  by  the  mismatch  between 
the  measured  field  and  the  MoM-simulated  field.  We  believe  it  is  this 
difference  that  drove  the  GA  to  zoom  onto  a  shape  that  is  similar  to  the  exact 
shape  but  has  a  lower  cost  value. 

Note  that  the  RMS  error  listed  above  is  particularly  high  for  the 
dihedral.  We  found  that  the  agreement  between  the  measured  data  and  the 
MoM-computed  data  for  this  shape  was  good  in  the  field  amplitude.  However, 
relatively  large  phase  difference  exists  (even  after  adjusting  for  the  rotation 
center)  between  the  two  results.  The  MoM  results  were  also  checked  against 
other  targets  in  the  Ipswich  data  set  and  the  phase  agreement  was  found  to  be 
good.  We  therefore  conclude  that  the  phase  data  for  Ips004  is  suspect.  As  an 


alternative,  we  also  performed  the  inversion  based  on  only  the  amplitude  of  the 
fields.  The  reconstructed  shape  is  shown  in  Fig.  3.  The  quality  of  the 
reconstruction  is  close  to  that  of  the  MoM-simulated  data  shown  in  Fig,  2(b). 

Conclusion:  A  genetic  algorithm  was  combined  with  a  moment  method  solver 
to  carry  out  shape  inversion  of  the  Ipswich  measurement  data.  We  used  a 
binary  bitmap  discretization  in  conjunction  with  a  geometrical  median  filter  to 
describe  the  shape  of  the  object.  It  was  found  that  the  algorithm  could  deal  with 
objects  containing  strong  multiple  scattering  effects. 
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Fig.  1  Ipswich  measurement  setup 


Fig.  2  Inversion  of  three  Ipswich  objects 

a  Real  shapes  of  Ips009,  Ips004  an  IpsOl  1 
b  Inversion  based  on  MoM-simulated  field 
c  Inversion  based  on  measured  complex  field 


Fig.  3  Inversion  of  the  dihedral  based  on  the  amplitude  of  the  measured  data 
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ABSTRACT 

Two  applications  of  the  adaptive  joint  time- frequency 
(AJTF)  algorithm  for  ISAR  image  formation  are 
presented.  First,  AJTF  is  utilized  for  ISAR  motion 
estimation  and  compensation.  Focused  images  from 
measured  radar  data  are  presented  to  illustrate  the 
effectiveness  of  the  algorithm  when  applied  to  in-flight 
aircraft  data.  Second,  the  AJTF  algorithm  is  extended  to 
detect  the  presence  of  chaotic,  three-dimensional  motions 
in  an  articulating  target.  Preliminary  test  results  on 
measured  data  show  that  the  algorithm  can  correctly  detect 
those  imaging  intervals  where  significant  three- 
dimensional  motions  exist. 

1.  INTRODUCTION 

High-resolution  inverse  synthetic  aperture  radar  (ISAR) 
imaging  is  a  promising  tool  for  non-cooperative  target 
identification  (NCTI).  The  main  challenge  in  ISAR-based 
NCTI  is  to  form  a  well-focused  image  of  an  articulating 
target  with  unknown  motion.  In  this  paper,  we  first  review 
the  application  of  joint  time-frequency  methods  for  ISAR 
image  formation.  By  using  an  adaptive  joint  time- 
frequency  (AJTF)  algorithm  to  estimate  the  phase  of  the 
prominent  scatterers,  we  show  that  the  target  motion  can  be 
estimated  and  a  focused  image  of  the  target  can  be 
constructed.  Results  of  applying  the  algorithm  to 
measured  ISAR  data  are  presented  and  discussed. 
Secondly,  we  report  on  our  recent  work  to  extend  the 
AJTF  algorithm  to  address  the  more  challenging  situation 
when  the  motion  of  the  target  is  not  limited  to  a  two- 
dimensional  plane.  In  particular,  we  discuss  our  research 
to  detect  the  presence  of  three-dimensional  motion  using 
the  AJTF  algorithm. 

2.  ISAR  MOTION  COMPENSATION 
USING  JOINT  TIME-FREQUENCY 
ALGORITHM 

We  first  review  the  application  of  joint  time-frequency 
methods  for  ISAR  image  formation.  To  form  a  focused 
image  from  raw  radar  data,  it  is  customary  to  first  carry  out 
a  coarse  alignment  of  the  data  in  the  range  dimension. 


followed  by  fine  motion  compensation  in  the  cross  range 
dimension.  Joint  time-frequency  techniques  have  been 
shown  to  be  a  useful  tool  to  carry  out  the  fine  motion 
compensation  [1,2].  We  assume  that  after  the  coarse  range 
alignment,  all  the  scatterers  are  located  in  their  respective 
range  cells.  The  radar  backscattered  signal  as  a  function  of 
dwell  time  r  in  a  particular  range  cell  can  be  written  as 

E{t)  =  y  4  expE-y— (/?(0 + cosOit) 

M  c  (1) 

+  y^sm0(t))] 

where  N  is  the  number  of  point  scatterers  in  that  range  cell, 
and  Ak,  yk  are  respectively  the  scattering  amplitude, 
down  range  position  and  cross  range  position  of  the  k*^ 
point  scatterer.  R(t)  is  the  residual  uncompensated 
translation  displacement  and  6(t)  is  the  rotational 
displacement.  Due  to  translation  and  rotational  motion,  the 
Doppler  frequency  versus  dwell  time  behavior  of  the  point 
scatterers  within  this  range  cell  is  not  constant  in  the  joint 
time-frequency  plane  (see  Fig.  1).  An  effective  JTF 
technique  to  extract  the  motion  parameters  is  based  on  a 
search  and  projection  procedure  to  represent  the  phase 
behavior  of  the  signal  E(t).  This  procedure  is  based  on  the 
adaptive  spectrogram  proposed  in  [3],  and  is  similar  in 
concept  to  a  one-term  matching  pursuit  algorithm  [4],  We 
shall  term  it  the  adaptive  JTF  (AJTF)  algorithm.  To  find 
the  motion  parameters,  basis  functions  in  the  form  of 

h(  t  )  =  exp[-j(  ajt  +  a2t^  )]  (2) 


are  chosen.  We  search  for  the  basis  function  over  the 
parameter  space  (a/,  a2,  a^)  that  best  represents  the  time- 
frequency  behavior  of  the  signal  by  maximizing  the 
projection  of  the  signal  onto  the  basis: 


max 
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E(t)h  (t)dt\ 


(3) 


After  the  time-varying  phase  for  the  strongest  point 
scatterer  is  found,  we  multiply  the  original  signal  by  the 
conjugate  of  this  phase  factor  to  compensate  for  the 
translation  motion.  This  algorithm  can  also  be  extended  to 
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Fig,  1 .  Fine  motion  compensation  is  carried  out  by 
the  Doppler  frequency  versus  dwell  time  behavior  of 
the  strong  point  scatterer  in  the  signal. 


Fig.  2.  ISAR  image  of  an  in-flight  aircraft 
obtained  after  AJTF  motion  compensation. 


multiple  range  cells  to  correct  for  higher-order  rotation 
motion.  After  applying  the  JTF  motion  compensation,  the 
standard  FFT  processing  in  the  dwell  time  domain  brings 
the  signal  into  the  cross  range  image  domain.  Fig.  2  shows 
an  exmaple  of  applying  the  AJTF  algorithm  to  measured 
ISAR  data  of  an  in-flight  aircraft.  The  shape  of  the  aircraft 
is  clearly  visible  in  the  resulting  image  after  the  AJTF 
motion  compensation. 

3.  THREE-DIMENSIONAL  MOTION 
DETECTION  USING  JOINT  TIME- 
FREQUENCY  ALGORITHM 

One  basic  assumption  of  standard  motion  compensation 
algorithms  is  that  the  target  only  undergoes  motion  in  a 
two-dimensional  plane  during  the  dwell  duration  needed  to 
form  an  image.  From  several  independent  examinations  of 
measured  ISAR  data  sets  recently,  it  was  reported  that  the 
presence  of  three-dimensional  motion  is  quite  detrimental 
to  focusing  the  image  [5-7].  We  shall  report  on  our  recent 
work  to  extend  the  AJTF  algorithm  to  address  the  more 
challenging  situation  when  the  motion  of  the  target  is  not 
limited  to  a  two-dimensional  plane.  In  particular,  we 
discuss  our  research  to  detect  the  presence  of  three- 
dimensional  motion  using  the  AJTF  algorithm. 

Allowing  for  arbitrary  three-dimensional  motion  in  space, 
we  consider  the  following  model  as  a  generalization  of  the 
model  for  two-dimensional  motion  in  (1): 

^  47Cf 

E(t)=j:Akexp[-j-^Xk+yk0  +  Zk0)]  (4) 

k=l  ^ 

where  0  is  the  azimuth  angle  of  the  target  with  respect  to 
the  radar,  and  ^  is  the  elevation  angle.  In  (4),  it  is  assumed 
that  the  translation  motion  has  been  removed  and  that  the 
standard  small-angle,  small  bandwidth  approximations 
apply.  This  model  reduces  to  the  standard  two-dimensional 
motion  model  when  ^and  ^are  linearly  related. 

In  general,  a  focused  image  cannot  be  obtained  from  the 
standard  two-dimensional  motion  compensation  algorithm 
when  three-dimensional  target  motion  is  present  due  to 
model  mismatch.  Therefore,  it  would  be  useful  to  detect 
the  presence  of  three-dimensional  motion  directly  from  the 
radar  data.  Our  approach  is  to  utilize  the  AJTF  algorithm 
to  extract  the  phase  behavior  of  the  radar  data  at  multiple 
range  cells.  We  first  parameterize  the  phase  of  the 
prominent  point  scatter  in  one  range  cell  using  AJTF. 
Next  we  repeat  the  same  procedure  at  another  range  cell. 
It  can  be  shown  that  when  the  target  undergoes  only  two- 
dimensional  motion  during  the  dwell  duration,  the  ratio 
between  the  parameters  («/,  ^2,  as)  extracted  from  one 
range  cell  and  those  corresponding  parameters  in  another 
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Fig.  3.  (a)  Simulated  2D  target  motion,  (b)  Phase 

behavior  of  the  prominent  point  scatterer  in  range  cell  1 
extracted  using  AJTF.  (c)  Phase  behavior  of  the  prominent 
point  scatterer  in  range  cell  2  extracted  using  AJTF.  (d) 
Ratios  of  the  extracted  phase  parameters  from  the  two 
range  cells.  Note  that  they  are  nearly  constant,  (e)-(h) 
Similar  to  (a)-(d),  except  that  3D  motion  is  assumed.  The 
resulting  ratios  in  (h)  are  no  longer  constant. 

range  cell  should  be  constant.  Therefore,  by  examining  the 
ratio  of  the  parameters,  we  can  distinguish  two- 
dimensional  motion  from  three-dimensional  motion.  Fig.  3 
illustrates  the  idea  using  simulated  point  scatterer  data. 
Figs.  3(a)-(d)  show  the  two-dimensional  motion  scenario 
and  Figs.  3(e)-(h)  show  the  three-dimensional  scenario.  It 
can  be  seen  from  the  results  in  Fig.  3(d)  that  the 
determined  ratios; 

Ci  =ai( range  cell  l)/ai( range  cell  2)  (5) 


(h)  C/=7.20 

C2=1.60 

C3=-L41 

are  nearly  constant  for  all  the  terms  in  case  of  two 
dimensional  motion,  as  expected.  For  three-dimensional 
motion,  the  ratios  are  not  the  same,  as  seen  in  Fig.  3(h). 

Fig.  4  shows  our  preliminary  results  of  applying  the  3D 
motion  detection  algorithm  to  real  radar  data.  Fig.  4(a) 
shows  the  degree  of  three-dimensional  motion  in  the  data 
for  20  different  image  frames,  detected  by  applying  our 
algorithm  to  the  raw  radar  data.  As  a  reference  for 
comparison,  Fig.  4(b)  shows  the  degree  of  three- 
dimensional  motion  for  the  same  20  frames  measured 
using  the  motion  data  derived  from  inertial  navigation 
instruments  carried  onboard  the  aircraft  during  data 
collection.  It  can  be  seen  that  our  algorithm  correctly 
detects  where  significant  three-dimensional  motions  exist. 
We  are  currently  fine  tuning  the  algorithm  to  achieve  faster 
and  more  robust  detection.  We  believe  this  detection 
algorithm  could  be  quite  useful  for  determining  the  “good” 
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(a)  Blind  detection  from  radar  data 
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that  the  algorithm  can  be  used  to  detect  those  imaging 
intervals  where  conventional  two-dimensional  motion 
assumption  would  fail.  We  are  working  to  devise 
algorithms  for  forming  focused  images  even  in  the 
presence  of  these  three-dimensional  motions. 
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Fig.  4,  Blind  detection  of  three-dimensional  motion 
from  real  radar  data,  (a)  Degree  of  three-dimensional 
motion  over  20  image  frames  detected  using  the 
proposed  algorithm,  (b)  Degree  of  three-dimensional 
motion  measured  from  on-board  instrument  data. 


imaging  intervals  from  which  focused  images  can  be  more 
readily  generated.  For  targets  that  exhibit  very  chaotic 
motions,  such  as  ships  on  the  ocean,  finding  such  intervals 
of  opportunity  may  be  very  critical  for  target  recognition. 

4.  SUMMARY 

In  this  paper,  we  presented  two  applications  of  the 
adaptive  joint  time-frequency  algorithm  for  ISAR  image 
formation.  In  the  first  application,  we  carry  out  fine 
motion  compensation  to  form  focused  ISAR  images  of 
articulating  targets.  The  AJTF  algorithm  is  used  to 
estimate  the  phase  of  the  prominent  point  scatterer  within  a 
range  cell.  The  higher-order  phase  error  due  to 
uncompensated  translation  and  rotational  errors  are  then 
removed  prior  to  the  image  formation.  Results  show  that 
well-focused  images  can  be  obtained  from  measured  data 
of  an  in-flight  aircraft.  In  the  second  application,  we  try  to 
detect  the  presence  of  three-dimensional  target  motion,  for 
which  a  well-defined  imaging  plane  does  not  exist.  A 
three-dimensional  motion  model  is  utilized  and  the 
linearity  of  the  phase  functions  of  the  prominent  point 
scatterers  between  different  range  cells  is  used  to 
distinguish  two-dimensional  from  three-dimensional 
motion.  The  AJTF  engine  is  again  used  to  extract  the 
phase  function  of  the  prominent  scatterer  within  each  range 
cell.  Preliminary  test  results  using  real  radar  data  indicate 


6.  REFERENCES 

[1]  V.  C.  Chen,  “Reconstruction  of  inverse  synthetic 
aperture  images  using  adaptive  time-frequency 
wavelet  transforms,”  SPIE  Proc.  on  Wavelet 
Application,  vol.  2491,  pp.  373-386,  1995. 

[2]  Y.  Wang,  H.  Ling  and  V.  C.  Chen,  “ISAR  motion 
compensation  via  adaptive  joint  time-frequency 
technique,”  IEEE  Trans.  Aerospace  Electron.  Syst., 
vol.  34,  pp.670-677,  Apr.  1998. 

[3]  S.  Qian  and  D.  Chen,  “Signal  representation  using 
adaptive  normalized  Gaussian  functions,”  Signal 
Processing,  vol.  36,  no.  1,  pp.  1-11,  Mar.  1994. 

[4]  S.  G.  Mallat  and  Z.  Zhang,  “Matching  pursuits  with 
time-frequency  dictionaries,”  IEEE  Trans.  Signal 
Processing,  vol.  41,  pp.  3397-3415,  Dec.  1993. 

[5]  V.C.  Chen  and  W.  J.  Miceli,  “Effect  of  roll,  pitch  and 
yaw  motions  on  ISAR  imaging,”  SPIE  Proc.  on  Radar 
Processing,  vol.  3810,  July  1999. 

[6]  A.  W.  Rihaczek  and  S.  J.  Hershkowitz,  “Choosing 
imaging  intervals  for  identification  of  small  ships,” 
SPIE  Proc.  on  Radar  Processing,  vol.  3810,  July 
1999. 

[7]  J.  Li,  Y.  Wang,  R.  Bhalla,  H.  Ling  and  V.  C.  Chen, 
“Comparison  of  high-resolution  ISAR  imagery  from 
measured  data  and  synthetic  signatures,”  SPIE  Proc. 
on  Radar  Processing,  vol,  3810,  July  1999. 


Conference  Paper  [11] 


ISAR  motion  detection  and  compensation  using  genetic  algorithms 

Junfei  Li  and  Hao  Ling 

Dept,  of  Electrical  and  Computer  Engineering, 

The  LFniv.  of  Texas  at  Austin,  Austin,  TX  78712*1084 


ABSTRACT 

Based  on  the  point  scatterer  model,  a  radar  signal  can  be  effectively  analyzed  using  the  joint  time-frequency  (JTF)  method. 
The  basis  functions  of  a  few  prominent  point  scatterers  are  believed  to  carry  target  motion  information  essential  to  the  ISAR 
imaging  process.  One  major  problem  with  the  JTF  method  is  the  computation  load  associated  with  the  exhaustive  search 
process  for  motion  parameters.  In  this  paper,  genetic  algorithms  (GA)  are  used  to  for  the  parameterization  process  in  the  JTF 
method.  Real  and  binary  coded  GA  are  investigated  and  their  performance  compared  with  the  exhaustive  search.  It  is  shown 
that  a  significant  amount  of  time  can  be  saved  while  achieving  almost  the  same  image  quality  by  using  real-coded  GA. 


Keywords:  inverse  synthetic  aperture  radar  (ISAR)  imaging,  motion  compensation,  joint  time-frequency  technique,  genetic 
algorithm 


1.  INTRODUCTION 

High-resolution  inverse  synthetic  aperture  radar  (ISAR)  imaging  is  a  promising  tool  for  non-cooperative  target  recognition.*’^ 
It  utilizes  the  motion  of  a  moving  target  to  generate  the  necessary  cross-range  resolution  in  forming  an  image.  Ideally,  the 
desired  target  motion  is  uniform  rotation  without  translational  motion,  under  which  a  two-dimensional  Fourier  transform 
converts  the  radar  data  in  the  (frequency)-(dwell  time)  domain  into  the  (range)-(Doppler  frequency)  domain.  Otherwise, 
motion  compensation  is  needed  as  an  intermediate  step  to  form  a  focused  ISAR  image.  Since  target  motion  can  always  be 
decomposed  into  translational  motion  and  rotational  motion,  a  typical  motion  compensation  algorithm  consists  of  two  steps. 
First,  a  point  on  the  target  is  focused  through  translational  motion  compensation.  When  there  is  non-uniform  rotational 
motion,  other  points  on  the  target  are  not  necessarily  focused.  Rotational  motion  compensation  is  then  applied  to  focus  these 
other  points. 

Recently,  we  proposed  an  ISAR  motion  compensation  algorithm  based  on  joint  time-frequency  (JTF)  analysis.^  The  basic 
idea  is  to  extract  the  phase  of  the  prominent  point  scatterer  within  a  range  cell  via  a  JTF  projection  technique.  The  data  is 
systematically  projected  onto  a  set  of  basis  functions  to  determine  the  one  that  best  resembles  the  joint  (dwell  time)-(Doppler 
frequency)  behavior  of  the  strongest  point  scatterer  in  the  range  cell.  We  use  polynomial  functions  as  the  phase  function  of 
the  basis.  Once  the  phase  functions  of  two  point  scatterers  are  found  by  JTF  analysis,  both  translational  and  rotational  motion 
compensation  can  be  achieved  with  good  accuracy.  However,  because  the  phase  parameters  are  obtained  by  an  exhaustive 
search,  the  algorithm  requires  heavy  computation  load  if  high  order  polynomials  are  used.  Therefore,  the  algorithm  is  not 
practical  in  dealing  with  cases  where  higher  order  target  motions  are  involved. 

In  this  paper,  the  genetic  algorithm  (GA)"*  is  explored  as  an  alternative  to  the  exhaustive  search  in  the  JTF  procedure.  GA  is  a 
stochastic  optimization  method  based  on  the  principle  of  natural  selection  and  evolution.  A  set  of  trial  solutions  forming  an 
initial  population  is  updated  iteratively  using  crossover  and  mutation  operators.  From  one  generation  to  the  next,  the  overall 
performance  of  the  population  is  improved.  Like  exhaustive  search,  GA  is  a  global  optimization  procedure.  However,  it  is 
capable  of  finding  the  optimum  in  a  very  high-dimensional  space  in  a  much  more  efficient  manner  than  exhaustive  search. 
This  paper  is  organized  as  follows.  Section  2  gives  a  summary  of  the  JTF  motion  compensation  technique  for  ISAR  imaging. 
The  concept  of  translational  and  rotational  motion  compensation  is  introduced  based  on  the  point  scatterer  model.  Genetic 
algorithm  is  described  in  Section  3.  The  general  GA  procedure  and  some  issues  specific  to  the  motion  parameter  estimation 
problem  are  discussed.  ISAR  imaging  results  using  some  real-world  data  by  the  proposed  JTF/GA  method  are  shown  in 
Section  4.  Conclusions  are  given  in  the  final  section. 


2.  JOINT  TIME-FREQUENCY  MOTION  COMPENSATION 


The  joint  time-frequency  projection  technique^  is  used  to  estimate  the  phase  of  a  prominent  point  scatterer.  We  begin  with 
coarsely  aligned  radar  range  profiles  over  a  certain  dwell  interval.  Within  a  fixed  range  cell,  the  data  can  be  written  as 


^  Ajrf 
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where  fc  is  the  radar  center  frequency,  to  is  the  dwell  time,  (p  is  the  target  aspect  angle,  and  r  describes  the  residual 
translational  motion.  The  point  scatterer  is  located  at  with  o;  as  its  complex  scatterering  strength.  Among  the  Nr  point 
scatterers  within  the  range  cell,  we  express  the  phase  behavior  of  the  strongest  one  as  a  polynomial  function: 


and  consider 


(2) 


/i(0  =  exp[-7P^(r^)]  (3) 

as  a  basis  for  the  radar  signal.  The  phase  parameters  are  then  found  by  searching  for  the  maximum  projection  from  the  radar 
signal  onto  the  basis  function: 

<  /i .  >=  arg  max  I  )dto  I  (4) 

In  the  search  procedure,  the  first  term//  can  be  obtained  using  the  fast  Fourier  transform,  while  all  higher  order  terms /2,/j,  ... 
are  obtained  using  exhaustive  search.  Figure  1  illustrates  the  process  of  JTF  phase  estimation.  Figure  1(a)  shows  the  coarsely 
aligned  range  profiles.  Figure  1(b)  shows  the  joint  (dwell  time)-(Doppler  frequency)  behavior  of  the  radar  signal  in  range  cell 
ro.  It  contains  three  point  scatterers.  The  tilted  trajectory  of  the  prominent  point  scatterer  1  implies  there  exist  higher-order 
terms  in  the  phase  function.  Figure  1(c)  shows  the  basis  function  h(tD).  During  the  search,  we  change  the  position  (/}),  tilting 
(/i)  and  curvature  (/>,  ...)  of  h  until  the  projection  of  h  onto  the  radar  signal  is  maximized. 

In  order  to  carry  out  both  translational  and  rotational  motion  compensation,  the  JTF  phase  estimation  is  needed  for  at  least 
two  prominent  point  scatterers  from  two  different  range  cells.  After  phase  estimation  of  the  first  point  scatterer,  we  multiply 
its  conjugate  phase  to  the  data  to  accomplish  translational  motion  compensation.  As  we  can  see  from  Figures  2(a)  and  2(b), 
the  trajectory  of  the  first  point  scatterer  in  the  JTF  plane  is  straightened  to  a  horizontal  line  after  the  translational  motion 
compensation.  This  point  becomes  focused  in  the  image  and  serves  as  the  origin  of  the  remaining  rotational  motion.  After 
phase  estimation  of  the  second  point  scatterer,  we  resample  the  phase  function  to  make  it  into  a  linear  function  of  dwell  time. 
As  shown  in  Figure  2(c),  the  JTT  behaviors  of  both  point  scatterers  become  horizontal  lines  in  the  JTF  plane  and  the  whole 
target  can  be  focused  in  the  image.  The  main  computational  issue  involved  in  this  algorithm  is  an  efficient  strategy  for  the 
phase  estimation  problem  in  equation  (4).  Exhaustive  search  is  very  computationally  intensive,  and  makes  motion 
compensation  of  targets  with  higher  order  motions  prohibitive. 


3.  GENETIC  ALGORITHM 

We  explore  the  use  of  genetic  algorithm  (GA)  to  improve  the  computational  efficiency  of  exhaustive  search.  GA  is  a  global 
optimization  method  based  on  the  concepts  of  natural  selection  and  e volution. The  stochastic  characteristic  of  GA  makes  it 
less  vulnerable  to  local  optimums  as  in  deterministic  optimization  algorithms.  Since  GA  is  based  on  a  set  of  population 
instead  of  one  single  possible  solution,  this  helps  the  algorithm  to  find  the  global  optimum  more  easily  for  a  complex  object 
function  with  many  local  optimums.  Another  advantage  of  GA  is  its  robustness. 


Figure  3  shows  the  flowchart  of  a  simple  genetic  algorithm.  It  consists  of  three  steps:  1)  Initialization.  Because  the  GA  result 
is  quite  independent  of  the  initial  state,  usually  the  initial  population  consists  of  random  individuals.  2)  Selection.  Selection  is 
based  on  the  fitness  of  the  individual.  Individuals  in  the  population  are  sorted  according  to  the  objective  function  evaluated  at 
the  corresponding  optimization  parameters.  According  to  some  stochastic  strategy,  two  or  more  parents  are  selected  to 
generate  children.  3)  Reproduction.  In  the  reproduction  stage,  crossover  and  mutation  operators  are  used  to  generate  new 
individuals  from  the  selected  individuals.  The  crossover  operation  produces  children  genetically  similar  to  parents.  This 
ensures  that  good  features  in  the  solution  are  kept  from  one  generation  to  the  next.  The  mutation  operation  generates  a  child 
that  is  different  from  the  parent.  Better  solutions  containing  features  not  included  in  the  parent  population  are  made  possible 
by  the  added  diversity  from  mutation. 

After  a  new  population  is  generated  by  a  series  of  selection,  crossover  and  mutation  operations,  the  selection  and 
reproduction  processes  are  iterated  until  a  preset  convergence  criterion  is  met.  There  are  several  choices  for  the  convergence 
criterion.  For  example,  if  we  know  the  true  optimal  objective  function,  an  RMS  error  might  be  set  to  ensure  the  solution  is 
accurate  enough.  In  case  we  do  not  know  the  truth  value,  we  might  set  a  fixed  generation  number  to  balance  the  solution 
performance  with  the  computational  cost.  A  third  option  is  to  stop  the  process  when  GA  has  reached  its  saturation  point. 
Typically,  the  rate  of  convergence  is  fast  in  the  beginning.  When  GA  has  run  for  a  long  time  and  the  solution  performance 
does  not  improve  much,  we  can  choose  to  stop  the  process. 

GA  operates  on  coded  parameters  instead  of  the  parameters  themselves.  A  coded  parameter  is  called  a  gene,  and  the  string  of 
genes  to  represent  a  set  of  optimization  parameters  is  called  a  chromosome.  While  in  traditional  GA  the  parameters  are 
usually  binary  coded,  real-coded  GA  has  also  been  developed  and  utilized.^  Binary  GA  utilizes  discrete  binary  bits  to 
represent  the  parameters.  Real  GA  utilizes  continuous  real  numbers  to  represent  the  parameters.  The  main  differences  in  the 
algorithms  of  binary  versus  real  GA  are  the  crossover  and  mutation  operators  depicted  in  Figure  4.  In  Figure  4(a)  depicting 
binary  GA,  each  chromosome  consists  of  8  bits.  To  make  the  crossover,  4  bits  from  the  first  parent  are  combined  with  4  bits 
from  the  second  parent  to  form  a  child.  To  make  the  mutation,  the  content  of  a  random  bit  is  flipped.  In  Figure  4(b)  depicting 
real  GA,  crossover  and  mutation  are  defined  as  simple  additions.  The  crossover  operator  is  a  linear  combination  of  the  two 
parents,  where  a  is  a  random  number  between  0  and  1.  The  mutation  for  real  GA  adds  some  randomness  to  an  individual 
chromosome  by  adding  a  random  vector  P  to  it.  In  case  the  mutated  result  is  out  of  the  search  space,  the  mutation  operator  is 
repeated  until  an  acceptable  value  is  found.  While  either  binary  or  real  GA  can  be  applied  to  an  optimization  problem,  it 
appears  to  us  that  in  problems  involving  continuous  parameters,  real-coded  GA  may  achieve  better  solution  performance  than 
binary  GA,  since  the  binary  search  space  is  limited  in  resolution.  Both  the  real-coded  and  binary-coded  GA  are  tested  for  the 
phase  estimation  problem. 


4.  RESULTS 

To  test  the  accuracy  of  the  GA  search  for  the  phase  estimation  problem,  a  phase  function  with  3^^  order  behavior  is  used  to 
generate  a  set  of  simulated  data.  Real-coded  GA  is  utilized  to  estimate  the  phase  coefficients.  In  our  GA,  we  use  50 
chromosomes  in  the  population.  The  crossover  probability  is  0.7  and  the  mutation  probability  is  0.3.  The  objective  function 
O  is  defined  as  the  projection  value  from  the  data  onto  the  basis.  Figure  5(a)  shows  the  convergence  curve  of  one  GA  run.  It 
is  seen  that  after  200  generations  the  projection  value  error  of  the  best  individual  is  3%,  while  after  400  generations  the  error 
is  about  1%.  Figure  5(b)  shows  the  truth  phase  function  and  the  phase  function  obtained  from  the  GA-searched  parameters 
after  400  generations.  The  GA  phase  result  is  nearly  indistinguishable  from  the  truth  phase. 

Next,  we  investigate  the  computational  complexity  of  the  GA  search  as  the  number  of  search  parameters  is  increased.  This  is 
shown  in  Figure  6,  where  the  computation  time  is  plotted  as  a  function  of  the  parameter  number  for  both  GA  and  exhaustive 
search.  As  the  dotted  curve  shows,  the  computation  time  for  exhaustive  search  increases  dramatically  with  the  number  of 
parameters.  The  solid  curve  is  for  the  real-coded  GA  search  and  it  grows  much  slower  than  the  exhaustive  search.  This 
means  we  can  use  GA  to  carry  out  the  search  much  more  efficiently  when  a  large  number  of  parameters  is  involved. 

Finally,  we  test  the  GA  search  for  the  JTF  phase  estimation  problem  with  a  real  radar  data  set.^*^  There  is  no  significant 
rotational  motion  of  the  target  in  this  data  so  the  phase  estimation  of  only  one  point  scatterer  is  needed.  Imaging  results  from 
real-coded  GA  are  compared  to  those  from  binary-coded  GA  and  exhaustive  search.  Figures  7(a)-(c)  show  images  generated 
for  a  typical  frame  using  exhaustive  search,  binary-coded  GA  and  real-coded  GA,  respectively.  The  corresponding 
computation  time  is  45  s,  22  s,  and  16  s.  In  the  20  frames  we  examined,  it  is  found  that  in  19  frames,  the  real-coded  GA  gave 
larger  projection  values  than  exhaustive  search.  The  resulting  image  quality  is  either  on  par  or  slightly  better  than  those 


obtained  from  the  exhaustive  search,  while  the  computation  time  is  about  1/3  to  1/9  of  the  exhaustive  search.  For  the  binary- 
coded  GA,  8  frames  are  found  to  be  of  inferior  image  quality  than  the  exhaustive  search  result.  Therefore,  real-coded  GA 
appears  to  be  the  better  choice  for  the  JTF  motion  compensation  problem  at  hand. 


5.  CONCLUSIONS 

One  problem  with  JTF-based  ISAR  imaging  is  the  time  complexity  associated  with  the  exhaustive  search  process.  In  this 
paper,  GA  is  used  to  replace  the  exhaustive  search  process  in  the  JTF  phase  estimation  for  higher  efficiency.  It  is  shown  that 
the  real-coded  GA  achieves  almost  the  same  image  quality  while  saves  computation  time  significantly.  With  GA,  the 
computation  time  complexity  for  higher  order  motion  is  much  smaller  than  that  of  exhaustive  search.  Therefore,  GA  is  a 
good  candidate  for  performing  ISAR  motion  compensation  of  targets  with  higher  order  motion. 
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Fig.l.  Joint  time-frequency  analysis  of  radar  signal,  (a)  Range  profiles  versus  dwell 
time,  (b)  Joint  (dwell-time)-(Doppler  frequency)  representation  of  the  radar 
signal  within  range  cell  ro  showing  several  point  scatterers.  (c)  Phase  of  the 
strongest  point  scatterer  is  found  by  projecting  the  radar  signal  onto  the  best 
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Fig.2.  Translational  and  rotational  motion  compensation  by  JTF.  (a)  JTF  display  of  the 
uncompensated  radar  signal  with  two  prominent  point  scatterers.  (b)  After 
translational  motion  compensation,  (c)  After  tranlsational  and  rotational  motion 
compensation. 
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Fig.3.  Flow  chart  of  a  simple  genetic  algorithm 
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Fig.4.  Examples  of  crossover  and  mutation  operations  in  binary  and  real-coded  GA. 
(a)  Binary  GA.  (b)  Real  GA. 


Fig.5.  Behavior  of  a  typical  GA  run  for  phase  estimation,  (a)  Convergence  curve,  (b)  GA 
phase  comparison  with  the  original  truth  phase. 
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Fig.6.  Computation  time  comparison,  (a)  Exhaustive  search.  (b)Real  GA. 
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Abstract 

With  the  assumption  of  rigid  body  motion,  traditional  high-resolution  ISAR 
imaging  is  capable  of  generating  a  (range)-(cross  range)  image  of  the  target.  When  the 
target  has  non-rigid  body  motion,  the  resulting  radar  image  does  not  preserve  the  spatial 
features  of  the  target.  In  this  paper,  we  investigate  the  extraction  of  micro-Doppler 
features  due  to  non-rigid  body  motions.  Joint  time-frequency  analysis  with  adaptive 
chirplet  basis  functions  is  employed  to  represent  the  range  compressed  radar  data.  After 
the  parameterization,  returns  from  different  moving  parts  on  the  target  can  be  easily 
separated.  The  algorithm  is  applied  to  analyze  the  measurement  data  from  a  walking  man. 
It  is  shown  that  the  rate  of  the  arm  swing  can  be  accurately  estimated  after  the  extraction. 

1.  Introduction 

In  inverse  synthetic  aperture  radar  (ISAR),  the  backscattered  data  from  a  moving 
target  is  collected  and  processed  with  a  motion  compensation  algorithm  to  generate  a 
(range)-(cross  range)  image  of  the  target.  Since  the  image  is  a  2-dimensional  projection 
of  the  target,  spatial  features  present  in  the  image  are  useful  for  target  identification 
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purpose  [1].  A  basic  assumption  in  ISAR  imaging  is  that  the  target  is  a  rigid  body  during 
the  data  acquisition  time. 

However,  non-rigid  body  motions  (associated  with  a  target  which  is  not  a  rigid 
body)  can  exist  in  many  real  targets.  For  example,  a  target  may  consist  of  multiple 
movable  components.  In  such  cases,  a  focused  ISAR  image  resembling  the  original  target 
cannot  be  obtained  by  the  tradition  ISAR  processing.  On  the  other  hand,  while  the 
geometrical  relationship  is  difficult  to  capture,  motion  features  are  still  embedded  in  the 
radar  data.  Recently,  the  term  micro-Doppler  has  been  coined  to  describe  the  small 
motion  features  of  a  non-rigid  target  [2].  Micro-Doppler  provides  an  additional  feature  to 
the  Doppler  return  from  the  main  body  motion,  and  has  been  studied  for  identifying  the 
natural  resonant  frequency  of  a  tractor-trailer  truck  [3]. 

In  this  paper,  we  investigate  the  problem  of  micro-Doppler  feature  extraction  from 
ISAR  data.  The  approach  we  use  is  adaptive  joint  time  frequency  (AJTF)  analysis  [4,5], 
with  chirplet  as  the  basis  function  to  represent  the  signal  [6].  This  work  is  an 
improvement  of  the  work  reported  in  [7],  where  chirp  basis  functions  with  constant 
amplitude  were  used  to  model  the  radar  return  from  a  helicopter.  We  apply  the  algorithm 
to  the  radar  data  from  a  walking  man  and  extract  micro-Doppler  features  due  to  the 
swinging  arms. 

2.  Adaptive  Chirplet  Representation 

Adaptive  joint  time-frequency  analysis  is  a  good  tool  to  analyze  the  time-varying 
Doppler  features  of  a  target.  Previously,  AJTF  was  utilized  to  achieve  ISAR  motion 
compensation  [8].  As  in  [7],  we  are  interested  here  in  parameterizing  the  radar  return  so 
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that  we  can  extract  the  different  components  of  the  signal  for  further  processing.  For  this 
purpose,  the  chirplet  basis  function  proposed  in  [6]  is  used.  The  radar  signal  as  a  function 
of  dwell  time  t  is  expanded  in  terms  of  basis  functions  as  follows: 

s{t)  =  exp[-(t  -  tk)^  lal  +  +  (1) 

;t=i 

where  c*  is  the  coefficient  with  the  k'’*  basis  function.  Each  basis  function  is  a  four- 
parameter  chirplet  with  time  location  4,  frequency  center /*;,  time  extent  Ofc,  and  chirp  rate 
fik.  The  parameters  for  all  the  individual  chirplets  are  found  adaptively.  The  chirplet  with 
the  highest  energy  is  extracted  first.  We  exhaustively  search  for  the  maximum  projection 
from  the  radar  signal  onto  the  set  of  basis  functions  in  the  parameter  space.  Once  the  best 
basis  is  identified,  we  then  subtract  its  contribution  from  the  radar  signal  and  continue  to 
search  for  the  next  best  basis.  This  process  is  iterated  until  we  have  a  good  representation 
the  original  signal. 

3.  Results  from  Walking  Man  Data 

A  particular  example  with  micro-Doppler  features  is  the  radar  data  of  a  walking 
man  [2].  Figure  1  shows  the  geometry  of  the  situation.  During  the  data  collection,  a  man 
walked  toward  the  radar.  At  least  two  different  motions  were  present  on  the  walking  man: 
the  steady  body  motion  and  the  swinging  arm  motion.  One  segment  of  the  radar  data  after 
range  compression  is  shown  in  Figure  2.  We  see  the  nearly  linear  range  variation  in  the 
figure  due  to  the  steady  walking  speed  of  the  man.  We  first  carry  out  a  coarse  range 
alignment  by  correlating  the  range  profiles.  The  signal  through  a  particular  range  cell  is 
then  used  for  JTF  processing.  Shown  in  Figure  3  is  the  spectrogram  (generated  using  the 
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short-time  Fourier  transform)  of  the  radar  signal  in  a  particular  range  cell  containing 
strong  micro-Doppler.  We  immediately  recognize  the  two  motion  components  present  in 
this  figure.  The  horizontal  trajectory  is  due  to  the  body  motion,  while  the  sinusoidal  curve 
is  due  to  the  arm  swing. 

We  apply  the  JTF  extraction  outlined  in  Section  2  and  extract  50  chirplets.  After  the 
parameterization,  we  separate  the  body  and  the  arm  returns.  As  can  be  seen  in  Figure  3, 
the  body  return  should  consist  of  chirplet  bases  with  both  small  Doppler  frequency /*  and 
small  chirp  rate  Pk-  By  applying  this  criterion  to  the  basis  functions  representing  the 
signal,  we  can  separate  the  body  return  from  the  arm  return.  The  results  are  shown  in 
Figures  4a  and  4b,  respectively.  Notice  that  the  main  features  in  the  two  returns  are 
preserved  after  the  JTF  extraction.  The  Doppler  features  due  to  the  arm  swing  shown  in 
Figure  4b  are  nearly  periodic.  This  period  can  be  easily  estimated  from  the  arm-only  data 
by  taking  the  autocorrelation  of  the  time  sequence.  The  result  is  shown  in  Figure  5a.  The 
peaks  in  the  autocorrelation  function  indicate  the  period  of  the  motion  and  it  is  estimated 
to  be  0.44  sec.  Note  that  without  using  JTF  analysis,  accurate  detection  of  the  swing 
period  would  be  more  difficult  due  to  the  presence  of  the  body  return.  The 
autocorrelation  calculated  from  the  unfiltered  data  is  shown  in  Fig.  4b.  As  we  can  see,  the 
peaks  are  much  less  prominent. 

4.  Conclusions 

Adaptive  joint  time-frequency  analysis  is  an  effective  technique  to  extract  micro- 
Doppler  features.  Using  the  measured  radar  data  of  a  walking  man,  we  showed  that  it  is 
possible  to  separate  the  return  of  the  body  from  that  of  the  swinging  arm.  Furthermore, 
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we  showed  that  the  arm  swing  period  could  be  accurately  estimated  from  the  arm  return. 

This  work  can  be  extended  for  other  targets  with  non-rigid  body  motions.  The  extracted 

micro-Doppler  features  may  be  exploited  for  target  identification  applications. 
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Fig.  1  Geometry  of  the  problem  —  radar  data  collection  of  a  walking  man. 
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Fig.  2  Radar  data  after  range  compression 
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Fig.  3  Spectrogram  of  the  radar  signal  showing  micro-Doppler 
features  of  the  swinging  arm. 


